{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elevation profiles extracted from SRTM over Mt Baker compared with ICESat-2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "import glob\n",
    "from topolib import gda_lib\n",
    "import topolib\n",
    "import gdal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read reference DEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dem_fn = '/home/jovyan/data/srtm_elevation/SRTM3/cache/srtm_wa_subset.vrt'\n",
    "# dem_fn = '/home/jovyan/data/baker_1_m/baker_2015_utm_m.vrt'\n",
    "# dem = rasterio.open(dem_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read ATL06 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = '/home/jovyan/data/nsidc/**/'\n",
    "ATL06_list = sorted(glob.glob(data_dir + \"*.h5\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract points from NED"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "ATL06_fn = ATL06_list[0]\n",
    "dataset_dict={'land_ice_segments':['h_li',\n",
    "                                   'delta_time',\n",
    "                                   'longitude',\n",
    "                                   'latitude'],\n",
    "              'land_ice_segments/ground_track':['x_atc']}\n",
    "\n",
    "ATL06_gdf = gda_lib.ATL06_2_gdf(ATL06_fn,dataset_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>h_li</th>\n",
       "      <th>delta_time</th>\n",
       "      <th>longitude</th>\n",
       "      <th>latitude</th>\n",
       "      <th>pair</th>\n",
       "      <th>beam</th>\n",
       "      <th>p_b</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>881.312744</td>\n",
       "      <td>2.522621e+07</td>\n",
       "      <td>-121.871270</td>\n",
       "      <td>48.700353</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0_0.0</td>\n",
       "      <td>POINT (-121.8712698055161 48.70035295193843)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>991.491089</td>\n",
       "      <td>2.522621e+07</td>\n",
       "      <td>-121.871498</td>\n",
       "      <td>48.701968</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0_0.0</td>\n",
       "      <td>POINT (-121.8714979934935 48.7019675199616)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1082.476440</td>\n",
       "      <td>2.522621e+07</td>\n",
       "      <td>-121.871683</td>\n",
       "      <td>48.703223</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0_0.0</td>\n",
       "      <td>POINT (-121.8716827998666 48.70322282685757)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NaN</td>\n",
       "      <td>2.522621e+07</td>\n",
       "      <td>-121.871883</td>\n",
       "      <td>48.704655</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0_0.0</td>\n",
       "      <td>POINT (-121.8718833403367 48.70465508248594)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NaN</td>\n",
       "      <td>2.522621e+07</td>\n",
       "      <td>-121.871913</td>\n",
       "      <td>48.704837</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0_0.0</td>\n",
       "      <td>POINT (-121.8719128466406 48.70483727673248)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          h_li    delta_time   longitude   latitude  pair  beam      p_b  \\\n",
       "0   881.312744  2.522621e+07 -121.871270  48.700353   1.0   0.0  1.0_0.0   \n",
       "1   991.491089  2.522621e+07 -121.871498  48.701968   1.0   0.0  1.0_0.0   \n",
       "2  1082.476440  2.522621e+07 -121.871683  48.703223   1.0   0.0  1.0_0.0   \n",
       "3          NaN  2.522621e+07 -121.871883  48.704655   1.0   0.0  1.0_0.0   \n",
       "4          NaN  2.522621e+07 -121.871913  48.704837   1.0   0.0  1.0_0.0   \n",
       "\n",
       "                                       geometry  \n",
       "0  POINT (-121.8712698055161 48.70035295193843)  \n",
       "1   POINT (-121.8714979934935 48.7019675199616)  \n",
       "2  POINT (-121.8716827998666 48.70322282685757)  \n",
       "3  POINT (-121.8718833403367 48.70465508248594)  \n",
       "4  POINT (-121.8719128466406 48.70483727673248)  "
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ATL06_gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pygeotools.lib import geolib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "ned_elevs = []\n",
    "\n",
    "lons = list(ATL06_gdf['longitude'].values)\n",
    "lats = list(ATL06_gdf['latitude'].values)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3860"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(lons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(ned_elevs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "ned_elevs = geolib.get_HAE(lons[0],lats[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(879.463)"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ned_elevs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ATL06_gdf['NED_elev'] = ned_elevs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ATL06_gdf['NED_ATL06_diff'] = ATL06_gdf['NED_elev'] - glas_gdf_aea_rgi['h_li']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ATL06_gdf['longtitude'], ATL06_gdf['latitude']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ATL06_gdf['longtitude']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df['geometry'] = glas_df['geometry'].apply(geolib.get_NED)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mgeolib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_NED\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlon\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m <no docstring>\n",
       "\u001b[0;31mFile:\u001b[0m      ~/pygeotools/pygeotools/lib/geolib.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "geolib.get_NED?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample DEM at ICESat-2 ATL06"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "ATL06_fn = ATL06_list[0]\n",
    "dataset_dict={'land_ice_segments':['h_li',\n",
    "                                   'delta_time',\n",
    "                                   'longitude',\n",
    "                                   'latitude'],\n",
    "              'land_ice_segments/ground_track':['x_atc']}\n",
    "\n",
    "ATL06_gdf = gda_lib.ATL06_2_gdf(ATL06_fn,dataset_dict)\n",
    "ATL06_gdf = ATL06_gdf.to_crs(dem.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly = gda_lib.dem2polygon(dem_fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi = gpd.sjoin(ATL06_gdf, poly, op='intersects', how='inner')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "points_xy = list(zip(glas_gdf_aea_rgi.geometry.x, glas_gdf_aea_rgi.geometry.y))\n",
    "\n",
    "refdem_sample = []\n",
    "for val in dem.sample(points_xy):\n",
    "    refdem_sample.append(val[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi['srtm_height'] = refdem_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi['diff'] = glas_gdf_aea_rgi['srtm_height'] - glas_gdf_aea_rgi['h_li']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD4CAYAAADrRI2NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2df5RU1ZXvv5umwRZhANOwsBEbSYsRQTr2EhgMKzHD4EjEEp8Kg1Ee82R0mbVC8DGhA++B78FAphPG5ZoVs0TzEp8tIoqlCSrhRZ288BoyTRpoCfYASqALBhgNymiHn/v9UbfwdvX9ce695/6q2p+1enXVubfqnqraZ9999tlnb2JmCIIgCOVDr7g7IAiCIESLKH5BEIQyQxS/IAhCmSGKXxAEocwQxS8IglBm9I67A2584Qtf4Nra2ri7IZQoO3bs+Hdmro76uiLXQpi4yXXiFX9tbS1aW1vj7oZQohDRH+K4rsi1ECZuci2uHkEQhDJDFL8gCEKZIYpfEAShzBDFLwiCUGaI4hcEQSgzEh/VI9izNNuOddsP4zwzKogwe8KVWJEZG3e3BCEQ2bYclr+2Bye7zgIABl1aiWW3j0GmvibmnpUOovhTytJsO57bduji8/PMF5+L8hfSSrYth0UbduHshc+zBv/xs7NYsH4nAIjy14Srq4eIriSit4loLxHtIaJvG+3riWin8XeQiHYa7bVE1GU69mPTe91IRO1EtJ+IniAiCu+jlTZmpa/SLvTk8OHDAHCNyHZyaNrc0U3pm2ncuDvi3pQuKhb/OQCPMvPviKg/gB1EtIWZ7y2cQEQ/BPCx6TUHmHm8xXs9CWA+gG0AXgdwK4A3fPe+TMm25eLuQknQu3dvAOhk5utEtpNB7mSX7bGusxci7Elp42rxM/NRZv6d8fgUgL0ALs63DMvmHgDrnN6HiIYBGMDMLZyv/vIsgEyAvpctj/18T9xdKAmGDRsGAJ8BIttJoZfMkyLBU1QPEdUCqAew3dT8FQDHmHmfqW0kEbUR0T8T0VeMthoAnaZzOmEaZEXXmU9ErUTUeuLECS9dLAv++NlZx+MyI/BOFLItcu2OjZdH0Iyy4ieiywC8DGABM39iOjQb3S2iowBGMHM9gIUAnieiAQCs7uWWPzMzP8XMDczcUF0def6s1LPklfa4u5AqopJtkWshKSgpfiKqRH5gNDPzRlN7bwAzAawvtDHzaWb+0Hi8A8ABANcgbwUNN73tcABHgn4AoSefnjkfdxfSBEFkOzXMWdsSdxdKApWoHgLwDIC9zLym6PBfAHiPmTtN51cTUYXx+GoAdQDeZ+ajAE4R0UTjPe8H8KqmzyEUIe4ed/LueFwFke1EoCKzWw98FEFPSh8Vi38ygG8CuMUUxnabcWwWei58TQGwm4h2AXgJwEPMXPi1HgbwNID9yFtLEvUQEstfkwVgN7Zu3QoAl0NkOxFI0EJ0uIZzMvNvYO3DBDPPtWh7Gfmps9X5rQCu99ZFwQ+FXY+CPTfffDMA7GDmhuJjItvR4xa0IOhDcvUIgiCUGaL4BUFIBAOrKpXOW5qVqLWgiOIXBCERjLmiv9J567YfDrknpY8o/hJGInuENPH/3leL2DnPsssrKKL4U4YXZS6RPUKaEH0eHaL4U4aXXbkS2SMIghWi+FOG7MoVBCEoovgFQYgdr+tREtkTDFH8JY4MECENeN21KwWHgiGKv8SRASKkAdm1Gy2i+FNGVaX8ZIIASLhyEESLpIxVM8fJjyYIkHDlIIgOSRmZ+hqsuXc8agZWgQDUDKzCpS6zALGMhFJEwpX9I4o/ZcxZ24IF63cid7ILDKD28ir8/cxxjq9p3Lg7ms4Jgk8qSIrtRoko/hQxZ21Lj0IUWw98hA2tzgu4XWcvhNktQQjMxKsHWbZXuNwPJGrNHyoVuK4koreJaC8R7SGibxvty4koZ1HAAkTUSET7iaiDiKaZ2m8konbj2BNGtSJBEbvqQ1KVSEg7rX84adnupiIkas0fKhb/OQCPMvOXAEwE8AgRXWcc+0dmHm/8vQ4AxrFZAMYAuBXAjwrl6gA8CWA+8iXr6ozjQgSIn78nhw8fBoBrxKiJn9PnrGel5y5IAp8wcFX8zHyUmX9nPD4FYC+AGoeX3AHgBaMw9QfIl6K7iYiGARjAzC2cL3b6LIBM4E9QJkxd847j8UGXOucylwiInvTu3RsAOsWoiY9sWw7jH/tl3N0oOzz5+ImoFkA9gO1G07eIaDcR/YSICk66GgDmhNmdRluN8bi43eo684molYhaT5w44aWLJcmElVuw7/injucsu32M43GJgOjJsGHDAOAzQIyaOMi25fCd9TsDy6b4+b2jrPiJ6DLk640uYOZPkLdwRgEYD+AogB8WTrV4OTu092xkfoqZG5i5obq6WrWLJcmctS04duqM63mZeid9JbgRhVEjBk13/u6lXdYKoAi32Wyz+Pk9o6T4iagSeaXfzMwbAYCZjzHzeWa+AGAtgJuM0zsBXGl6+XAAR4z24Rbtgg1WUTxWFErW9e0t8fx+iMqoEYPmc+asbcGZ82r+e7fZrKwCeEclqocAPANgLzOvMbUPM512J4B3jcevAZhFRH2JaCTy/s7fMvNRAKeIaKLxnvcDeFXT5yg5lmbblaN1ls/ID4zv3yXx/D4giFETKaoGDZA3amQ2qx8Vi38ygG8CuKUoyuEfjCiG3QC+BuA7AMDMewC8COD3AN4E8AgzF5LIPwzgaeR9owcAvKH105QQXsLUCgPDbYBIPH938u54XAUxaiLFS/hxwajp16fC8TyZzXqjt9sJzPwbWE9lX3d4zUoAKy3aWwFc76WDgjP3TRwRdxdSy9atWwHgchhGjdH8PQCziWg88l6EgwD+FsgbNURUMGrOoadR81MAVcgbNGLUaKBgzKy8cywWrN9pe96jL+6UmYEHXBW/kGxWZMZ2e96LAKfQ56XZ9h6vKVduvvlmANjBzA1Fh8SoCQm/lnmmvsZR8SsuFwgGkrIhgcxZ2+L7tX89wXkGIDsdhTjxUjPazb1TjLh71BHFnzCybTllH+jkUYN7tIk1LyQZLzWjV97ZXZZlk6I+RPEnjO++rBZ5UzekH5ofnBRybwRBHxNWblE+d/KowT189rJJUR+i+BOGXc4SM4MurcT+459i8uq3LKe3suArJI1sW05pIyKQl18ro0YWb/Uhi7sppFCfNHeyCwtNC15Oi1+CECdNmzuUz31u2yFs2fNv2L5kKgBvcf+CGqL4U84FAIs27ISXEP1sW06sJyFScie7PJ1/7NQZTFi5BV8ccpko/RAQV08J4HVflhfrSxB04CdJ9bFTZzwrfYnsUUMUf4KIKsugV+tLEILCEcXZS2SPGqL4E0SUWQbFMhJKEYnsUUMUf0KYuuadSLMMPvqiLAQL0eBWREiIHlH8CWDO2hbXQiu6kS3uQhSoFBHSjcxm3RHFnwDiilqQykVCmHiJ3dfJog0ym3VDFH/MxKl8JW+PECaqu9B1I9nH3RHFHzPPbxflK5QmKrvQw0LcPc6I4o8ZpxTKgpBW4la84u5xRqX04pVE9DYR7SWiPUT0baO9iYjeMwpSv0JEA432WiLqMlXr+rHpvW40qnbtJ6InjGpFZUvcgwMIlgI67Rw+fBgArhHZ1o+X9MthIO4eZ1Qs/nMAHmXmLwGYCOARIroOwBYA1zPzOAD/CqDR9JoDzDze+HvI1P4kgPnIl6yrA3Crjg+RVhYmIKSynLfD9+7dGwA6Rbb14yX9clhI8II9roqfmY8y8++Mx6cA7AVQw8y/ZOZzxmnb0L3YdA+MOqYDmLmF88VOnwWQCdT7FJNtyyXGzZOEmUccDBs2DAA+A0S2dZIUhSvBC/Z4StJGRLUA6gFsLzo0D8B60/ORRNQG4BMAS5n5/wKoAdBpOqfTaLO6znzkrSeMGFGaKYZ1TYUPrp7eo+2LjZtwzsNNZdEGqVcahWyXg1wDehRuLwLeX9VdtiVLpz6UF3eJ6DIALwNYwMyfmNqXIO8OajaajgIYwcz1ABYCeJ6IBsC6YLulemLmp5i5gZkbqqurVbuYKsKcCu9fNR0D+qqXrSt3f2hUsl0Ocq1r9mg1G25+cJJl1TknkjL7SBpKip+IKpEfGM3MvNHU/gCAbwCYY0xxwcynmflD4/EOAAcAXIO8FWSeMg8HcETHh0gbUWxh3/3YrTi4erryQClXdw/ySltkWxP/dcMuLe9TYbM23vzgJBxcPV252JC4e6xRieohAM8A2MvMa0zttwL4LoAZzPyZqb2aiCqMx1cjv9D1PjMfBXCKiCYa73k/gFe1fpqUoGsLe92Qfq7nFAaKGwvLsIiLoc+vgsi2Ns5pWriaPeFKx+MrMmNxcPV0pTEg9ETF4p8M4JsAbjGFsd0G4J8A9AewpSi0bQqA3US0C8BLAB5i5oJj7mEATwPYj7y19IbGz1J2bFn4VeVz3Vw/5ejt2bp1KwBcDpFtLehyq0weNRgrMmPdT4TaGBB3T09cF3eZ+Tew9mG+bnP+y8hPna2OtQK43ksHSw1dLpXH7x3v6fzdj92K2sWbtFy7VLj55psBYAczNxQdEtn2QbOGXeiVvWBZbzcIz207pHwjKRdk527EPPZzPYUiwojCKWM/v6ABHcVWmu72ZtAAQFWlqDGvyDcWMYVC6UFQXdgqpl8fZ3ePVC8S4saPQbNq5jjXc8So6Y4o/gjRJXx+p60r73R+nVQvEvyiQ7a9ui8LqNwsxKjpjij+CNHh5rELc1NBZYCIZST4QUcK5jA3EYpR0x1R/BGiw83jFuYWFLGMBD8ETcE8tH+fQK/36/4sV0Txp4i6If1Cj04Qy0jwio5Z4vYlUwO9XqJ2vCGKPwUQ8haNl7h9QYiKIC7Mof37KG0w1MGElVsiuU4a8JSkTfDO0mw71m0/jPM+Yt0IwAcRDQozE1ZuCWyBCaVNti2Hps0dOHKyyzrhlgtRKXszx06dQbYtV/YJCQGx+ENlztoWPLftkC+lDwBzQvBb1gyscj2nMEAEwYpsWw4LX9yJnE+lHyeNG+OpA5w0xOIPiWxbLnAK2TD8loumjcYChbw8jRt3h24ZZdtyln25b+II8dkmmO9t3B2olsTAqkp9nfFIV0SpaJdm23skiKsb0i8x7lqx+EMiaL59r+lnVcnU16C3QkRo2APETukD+S32kl8luXwWUDaWzxijqSfdUR0zYcuWldIH8skZxy17M9RrqyKKPySC5tvXna/EzP5Vav7VMN0933GZdUg63dKkbki/0GaSqmMmbNlyev9PTp9PhFEjij+B+N3BqJuw/KFLs+2p8w0Legjb1RF3PP+ctS2u5+hIZhcUUfwJI0yLyEzvXu7+nrDcPc8nQPAFfwSZBUZh0KiuDYU1m1VZ19ORzC4oovgTRlSLPz+4+wal88IYIEkpMi94x+/a1eRRgyMLo1Swacqy8JAZlQpcVxLR20S0l4j2ENG3jfbBRLSFiPYZ/weZXtNIRPuJqIOIppnabySiduPYE0a1opLDr7KMMtpBdRCWcvjb4cOHAeAakW11/K5dhblmVcyae9xnFheg36hJUwi0isV/DsCjzPwlABMBPEJE1wFYDOBXzFwH4FfGcxjHZgEYA+BWAD8qlKsD8CSA+ciXrKszjpccfq2isKId7HBL0wxEF/4WB7179waATpHt0kLVqNFt9euqtREFroqfmY8y8++Mx6cA7AVQA+AOAD8zTvsZgIzx+A4ALxiFqT9AvhTdTUQ0DMAAZm4xilc/a3pNSeHXKop6R6FbmuYCdpbMnLUtqF286eKfysKWF6sobAtq2LBhAPAZILJdjjiZNEuz7RjV+DpqF2/CqMbXlSJxvCRhjHt24MnHT0S1AOoBbAcw1CgyDeP/EOO0GgCHTS/rNNpqjMfF7VbXmU9ErUTUeuLECS9dFDygeqNZtKGnZTRnbUuPhaytBz5yVf5eZkPfi9DNFIVsi1wnDysFXLzj/jyz9r0lcWfBVVb8RHQZ8vVGFzDzJ06nWrSxQ3vPRuanmLmBmRuqq6tVuyiEhJW3xy56wS2qwctsKOhGIVWiku1ylWuVNCG6UXFjAj0VsNOOe53RaHFnwVVS/ERUifzAaGbmjUbzMWOKC+P/caO9E4A5afxwAEeM9uEW7QLi3cbulbinqZohiGyHyqJpoyO/pqobs1gBO81ISykaTSWqhwA8A2AvM68xHXoNwAPG4wcAvGpqn0VEfYloJPILXb81psyniGii8Z73m15T1vSi6Bd2C/jZ8OK2iJWWG0PeHY+rILKthF9XRxzZMP1eM+iO+7SgYvFPBvBNALcQ0U7j7zYAqwFMJaJ9AKYaz8HMewC8COD3AN4E8AgzF77NhwE8jfyi2AEAb+j8MEnAj9KroHgGB+Bvw4vbIpZO/6XKgrFftm7dCgCXQ2RbiWZJo6GVOA0k1+yczPwbWPswAeDrNq9ZCWClRXsrgOu9dDBt+FmQTEPEpJdsnTr9l0EznDpx8803A8AOZm6wOCyyXUTaPB0E/X22y+fvR4lHkQHXDtm5q5moFiSjppTj+YXSRLWeRWFWqaK87WazfvbuxDmmRPELqBvST+m8JGQVFARVVN2YhVmliovSbjabtrUBUfyCcn6g57YdwtQ17/i+jtw4hKRSu3hT7CGWUSKKX/DEvuOf+n7tuu2H3U+yIC1RQqWM399g0KXxhinHnaY5qYjiF7RjN9j91h5OUw6UUsXvLurp44Zp7ok3SrGEZ7Yth8mr38LIxZswefVbvm7KUnNXI2l2ZVT20hNdVFlBWHa79z0J/fpU2PpJveRAKVCcUqJv7174/l3jYouiSDt+gxZ+setoSSrfuMi25fCd9TsvRivlTnZdrGbnRbbF4tdImssFNt2tp0hG03+6wZdyVd1pqYJVHqHT5y5gwfqd4jaKmHLwm0cpU4++uLNHiCoD+LuXdnl6H1H8mgjy46vmFQkTXZZw0+YOX99Fpr4GVZXW4ug1nYVT7H/cybHSSJAF/SQQ9jpD0+YOx+MVNqUZ7NrtyLblcN7GW3rG7oANovg18d2X/WeS7EpZKJgTuZNdaNzY7kv5r5o5rodA9oLedBblYIHqJsiCfhLw43r0Qu5kl+Px2ROu9NRuh85UKaL4NXH6nH8HeRjVgPyga+bRdfa8qxVkR/G36PVbTcL3WEqUwvepc13HykZ3M9wbrhpsWQ7ygxP/4enabmtdXsacKP6EEGXueTt0+tmLrSCnnDuFkDs7P6UX/6XfG45gjXyf3bFyqDgFqw2sqsRjP99jmdlz64GPtAaEHHGZeZgRxZ8QkpDqIcyIFye/e8NVgwHY+ylV/JeFEDe3abfgDR3fZxJmDTrX0cyfx+2zLZ8xxtFSVwkIKVQDc+MKD3UPyjKc0yrq476JI9Bw1WA0be7AkZNduGJgFRZNG62kDNMcxhkFbt9P0+aOQDedbFsOj27YhfOllDDdJxNWbsGxU2cuPh/avw+2L5kaY4/ytW3jDqNdeedYLNBUY3f5a3sufh63HD2Z+ppA17XSVXZ4qXtQdorf7ot8btuhbnff3MkuLDJcDG5CqyuM0y7zX5QMurTSV9y8E26Vi7xMUa1Y8kq7KH0A45a9iU9Odw8UOHbqDGoXb0IFEc4zo4IIsydcGWlsfWENK07ZDqqAzZgDBFRy9AysqvQVVOBUDcwKieN3wMsXefY8R7pr1GssbhiEEQHhppNVpqhOswYvCbLiTiEQFtm2XA+lb8ZP/VidtRCSsPu6wlv0pDb8RqV5+c769vamylUqcP2EiI4T0bumtvWmwhUHiWin0V5LRF2mYz82veZGImonov1E9IRRqSjxuFm/45a9qe1aXmNxw0CnVabq21WZouqYVfWi8EP74sKrYnX7Pr1am27onkX64Yf36NmkWEDVxet3THn5zr5/1zhP761ym/gpgFvNDcx8LzOPZ+bxyNcr3Wg6fKBwjJkfMrU/CWA+8uXq6orfM604WVl+SMJ6ga7EVqqRSoWBEbYlsOae8d0G4bx58wDghlIwavwoVqcbcxIsdN3odjX5TTpYzISVWwK/h9fP5qr4mfnXACxv/YaA3wNgndN7GAWrBzBzC+cLnT4LIOOppzHhtGv02iXuK+1eSUJ5O13+38/OXnC1+mtMbh63whlWN0UvN8riwTF37lwA2GduKyejxilUMwkWehhMHjVY23v5TTpYjHkxPiqC+vi/AuAYM5sHz0giaiOifyairxhtNQA6Ted0Gm2Jx8k/96cQXDPxO3vy9LbaceIDtwU1s5vH7Ybz3LZDqF28CbWLN130P7stHDsxZcoUADhndSxtRs2lNukunLAL1UzCrDMsmh+cpOV9vIaouo2mglyPXvqG5/f2mtIECK74Z6P7wDgKYAQz1wNYCOB5IhoA689tq+OIaD4RtRJR64kTJwJ2sTuSn1uNH9x9QyTX8Tv93nrgI4xb9qbrwnEBH7+7dqMmTLnu09t7nLqdMgor2aDOxeK4UYkQMt+MVctA+kkm6Gfx2LfiJ6LeAGYCWF9oY+bTzPyh8XgHgAMArkF+MAw3vXw4gCN2783MTzFzAzM3VFdX++2iJSsyYz2t7ttNh3X45ZJM3GGlKnhZX/HhvtJu1IQp1x/7CBeMenapc7E4CDrdPU78/czPF1y9yp+X0FM/YzWIxf8XAN5j5ovWDhFVE1GF8fhq5P2d7zPzUQCniGiiMYW+H8CrAa7tm2xbDr08uDGspsNfbNwUi1+u1LDKxqnJwxSIMI2asPCya9OObFsOtYs3aehNstHl7nFi8qjBkRhPfoeLSjjnOgAtAEYTUScR/Y1xaBZ6+j+nANhNRLsAvATgIWYu3OYfBvA0gP3ID5o3fPY5EE2bO3DWo2/ePO0at+xNnAtoKrn9WEnY4g6EH2WzambPELS/nqDfFefDwkudUfO1a4PNILJtOS0bnKKyppOO1c3Fa6y9CqoupGJUonpmM/MwZq5k5uHM/IzRPpeZf1x07svMPIaZb2DmLzPzz03HWpn5emYexczfMhbCIsfPLtFCDne3TTKquP1YjQlI2AY4uwJ0rJVYWUS6d5ROHjXY1sKbPXs2AFyLEjBq3n4v2JqBjvDNuiH9XK3ppBg1Tjx+r954/wJeY+3duG/iCN/jpexSNlwxsMpz4qnCdmsdg6PwY23afdQ2ZK4rAQnbgHyopd13FXQB0CkIRUfaiKrKXtj7P//K8Zx169bhhRde2M3MDeZ2Zp5bfC4zv4x8eGcPmLkVwPX+exscv8nUCqkUgn7fA/pWYMvCr7qe17hxd+LXj4JmyrWb9ehKG3Fw9fTA71F2KRu+dm11DxdGVaV7RMTVizcFHhwHV0+/eId220GaBMvIS9InrziVLzipIYbcyo1UyvjdMrZg/U4tfv3dj6ltXUiKUeNE0Ey5UawhBKWsFH+2LYf1/3K4hwvjrhtrXBcVdYurm9WzaIOehFJBCNMyc1qM1LFQmXSrUjfxOE6FqKkb0k/L+5SV4n/s53ssF3Y37T4ayqJiEFJgGPmmohc5ziaCzjRkr0ayKaV4/mKG9u/jeDxo1JqKO02FslL8dq6aP352NtI0taqkvci1HT+8+wZHizxTX+M6gOzoTfoXiNOOziIkKrhlQN164KPYXZl+druq0HjbdY7HgxiYuqx9oMwUf9rYd/zT2AdIGCxYv9M1LcD2JVNxicc8ugP6VmD/quALX6WGzpKaVhTPsFQyoMadgtxvqmQ33HbdrsiM9TUjrRvST5u1D5RhVI8T/fpUeMrtHgVBq1MllUJUkJN1/t7K2ywL5ziFaAo90VmExIri31DlenGnIA/zO1ngUnFsRWYsVmTGYuqad7Dv+Kfdjg26tBLLbh8T+pgXxW9CZ3m2Yqymliphi3HXkL1v4ojQcrc8t+2Qq1tGFLyQRlQqjum04L0irh4TYd5lraaWKlPiuDMYiL88+aQ1m2Za+61CUjZh2iGKPyKsbioqN5pSj9Ir5QiPqHAqCBJXCgWV6JXnth2KdQ3LRyZrZZK+X6GsFL/dSn5YK/yCO0nJ2JhmnAqCFFxlYdUattPvqtErj74Y336VprvDSc1QIMlGTVkp/uUzxqCyyBSp7EWhrfCroBJql2QB0kEpRi5FSYXNtl1ze1i1hu3yTqm6CONc4w17ATXJRk1ZKf5MfQ2a7r4BNQOrQMjnomkqiikPI+bZabqtEmoXtwCFvSGqkARP8MfsCVe6toeh5C6pIFkDciGpRk1ZKX4gPwC2Lr4FH6yejq2Lb+kxIHTHPLuFHmbqaxKfyjbswX3SRxER4XMKseEFC7+CyDJzo860wEP798F7K29zPCcs95JOwu5jUovWSzhnEZn6GrT+4SPfIYx+YsybH5xUFgUwhPAoxIY78f27xvkOV+5FwJp7xnuaOSy7fUyo+wd0EHYfk1q0vuwsfhVWZMb63h5d7nHnhHw+8zAjJgR/ZOprfOeaf3/VdM/uItXZbJxhnV4+0+P3jk/87FwVlQpcPyGi40T0rqltORHliGin8Xeb6VgjEe0nog4immZqv5GI2o1jTxjVikIh25bD5NVvYeTiTZi8+i1ffrb9RTvqVAgybXSbhsftK1R1E8yZOAKZ+hrs+/vpOLg6//f4veMTOe2fN28eANxQTrId9S5wFUPIKRw1SWTqa9D84KRucl2jIZNsHKiM5p8CsEq2/Y/MPN74ex0AiOg65KsXjTFe86NCuToATwKYj3zJujqb9wxMti2Hxo3tyJ3sAiO/87VxY7vnAeIn2CBI5IRbdZ6gxSGColo9yMrdkKmvQdt//0vb18S1SW3u3LkAsM/iUGJle9GGXd1ke9GGXZ5kOw7r2m2W4RSOGgV+AzoK64VJNGrcUCm9+GsAqmEldwB4wShM/QHypehuIqJhAAYwc4tRcvFZABm/nXaiaXMHus52z7fTdfY8mjZ3KL+HX+s6iDXl9tqgxSGCosNStJsm/3lM0+cpU6YAwDnF02OX7eWv7cHZC92V5NkL7CkqKqz0G04kPdeUSkCHk3ESVqhsmATxxH6LiHYbrqBBRlsNAPO8rdNoqzEeF7dbQkTziaiViFpPnPBWS9Supq5qrV2/Rad1WK1Rp8+Nmj1HTnlqj5FQZDuIXGaKEZ0AABcJSURBVAP20U8nu84quYDGLXvT8zXLgaA76JN+Y7PCr+J/EsAoAOMBHAXwQ6PdSv+xQ7slzPwUMzcwc0N1dbWnjtlVb1Kp6uRX6QN6UiuEnT43bpwUV4IITbaDyLUbbi6gqWvewSen/WWe1bGPI62+8AJp738xvhQ/Mx9j5vPMfAHAWgA3GYc6AZh3kwwHcMRoH27Rrp1F00ajsiiPe2WFc8WnAgsDbB+32z3pBSfLIe5kbYB7H4L0Me7F6wJJlm0nX7KTCyjbluuR/tcLOvZxhFm/OQqC9D8psm3Gl+I3/JoF7gRQiIp4DcAsIupLRCORX+j6LTMfBXCKiCYaEQ/3A3g1QL+dKba3TM/tpsRLs+24EMBsD3uBKgnJ2uy25xdw66OT4krK7t0ky/ay28f0SH7mlAytMJNa8kr8WTDT6A4xE6T/Sfj+i3HdwEVE6wB8FcAXiKgTwDIAXyWi8ciP9YMA/hYAmHkPEb0I4PfIL5o9wsyF+eXDyEcIVQF4w/jTTtPmDkvrp2lzBza0HuqW/qAQ8QMAzQEXvXRNBWsGVlnm4E9C5MCKzFjHxUG378Bps0wc7p7Zs2cDwLUAKA2y3fqHj3oYJyrGStKKCyWRyaMGh5Ya5dMz55Xy80eJSlTPbGYexsyVzDycmZ9h5m8y81hmHsfMMwyrp3D+SmYexcyjmfkNU3srM19vHPuWEQGhHbtF3NzJLssfthDxE7QzuqayX7vW2vf7sbGAl2TcvoMkCT4ArFu3DgB2p0W244jIiYq4ZTvsjZdeogqjoOT2V6os4hajGvHjhC6l9vZ71tEeFzgZwmOXwrqqspfSd+C0DhD34E8yfjO0pqXYSRJcfSpZTv0SdyW9YkpO8dtZzE5ckqD8Ak43oSQIj11q61Uz1TZ4OZnCC190LlRdzvh1Qzy/PdgsoSqisZGEyC6VLKdBSJJsJ0fjaSDblvM1Hf5TwM1ROseGnxlLlKiktvZLUmY1ScEciOCXIAELAJRv6KWAapZTK1TW4JKUqbOksnP6URqVvYCgm2J1VvJZNG104jMaZuprQvPXJ2FWkwQKqUeKd6F7IeiCZd2QfolblwkblSynVqhk+UxSps6Ssvj9+OqDKm3dg6PcBpoVfhPrlRJWqUeiZsvCr2p9P5U6vGklbeO2pBS/HzfJ4pd2Bbqm7sFR6qgMfr+J9UoJHTOfINZ+GFXXgrqdSoGkyHRJKX4/IZV/ClD0029uczfsUiDrrKAUF6pFuL0m1is1dESSBCGMqmtO+zzCLu8ZBSo1PJJi0KRek8xZ24LaxZtQu3gTFqzfiaH9+0R27bCmd9+/a5zlDk3V1MhJZkVmrHJqBx1htmlkabY91lTFYeWlWTRtdI+IsFJCZfafFIMm1Yp/ztqWHtPZY6fOYGj/Pt2iTtJWNSdTX4M193QvXjLgkvh37upCVaUlPcIpDJZm210j0wpyHdaMM6y8Opn6mh55tAo8t+1QIizhKEhCAENqFX+2LWfrwzx26gwWTRt9saB6y/v+fJ2P3zveNo45ivjmj01RACe7zmLh+tKIc1dxY1RVVqQ+sZcfVFKHFDJwhhX9FeZCpVNdibiLDelAJawzbjcekFLFn23LYaGL0C9/bc/FGGi/i0qZ+hrcdeNwy2N27bpo3LgbxUPkgtGedtw2xBABq2aOTV2khA6icPDY7b6Om7iLDelApShL3BXHgJQq/qbNHT2UYjEnu85ezE/uh7oh/ZBty9laYHapFXTRZTMI7NrThNvCIXP6wuPSxDduGOZ+kuALFblNQm7/VCr+KBb9HvlaHRau32l74yjXhUddlEIUR1rZuKPT/SQhNJLgwkyl4o9i0c9tVlGOC486CSNcMO1EkVBtYFVlrC6VJKQXDxu3z5iE2WwqFf+iaaND77ibRZ+Eu7ZQOji5FXWyfEa8hcHTWJjcK2n4jK760yg4fZyI3jW1NRHRe0ZB6leIaKDRXktEXUS00/j7sek1NxJROxHtJ6InjGpFvsjU12CUwmYJv9w3cYSrRR/2Xdsu3LmEw6AjZ968eQBwQxJkW0dNCBXc5FZlE1KY1y8F0vAZVQznnwK4tahtC4DrmXkcgH8F0Gg6doCZxxt/D5nanwQwH/mSdXUW76nM0mx7oBqibqzIjI3dorfb4aq68zXtRBG2OnfuXADYV9Qci2xHFdvt9r1GkYLEbnEzCYue5YJKBa5fA/ioqO2XzHzOeLoN3YtN98CoYzqAmVuM6kTPAsj463JpVyIqECRFbFpw2goRRQrbKVOmAPkyiheJW7bD5L6JIxKxa3TRtNGoqqzo1lZO+zaSUBxHh6t8HrrXGB1JRG1E9M9E9BWjrQaAOZSg02jzTNiWYCHaJAkDZEVmLA6sug0HV0/HgVW3lZTSB5wzoyYkhW2ksh02KzJjExGNlqmvwaqZY7vtri+nfRtJMFwD5eMnoiXIW0zNRtNRACOY+UMiuhFAlojGwLrinq1Lk4jmIz91xogR3V0bYZdoKyhXpwGS1A0waSNTX5PY2gNhyLaTXIdNoZN/VlVpW+0qyuWjMGs6JIGagVWJSM1gh2+Ln4geAPANAHMKxaWZ+TQzf2g83gHgAIBrkLeCzFPm4QCO2L03Mz/FzA3M3FBd3b2UYpgl2syC71SOMe7IiFLiUpvv2a49CsKSbSe57tenwuol2phjzGTPnrcP5Yx/P2np4KcEbJT4Gl1EdCuA7wKYwcyfmdqriajCeHw18gtd7zPzUQCniGiiEfFwP4BXA/deM3NMm4qcyjGWsqUSNXaKyElBhUlcsr3yznDdeIWZ7Kdn4i3uUi6EvbM/KCrhnOsAtAAYTUSdRPQ3AP4JQH8AW4pC26YA2E1EuwC8BOAhZi4sDD8M4GkA+5G3lsy+00Rg9qGL9RMNdvfXKPYYzZ49GwCuRQJkO0xjQjU7rbgw9ZGEtRQnXH38zDzbovkZm3NfBvCyzbFWANd76l2ESAqB8mPdunV44YUXdjNzg6m55GS7+cFJSueJC1MfVzj4+JNwg03lzt0wMFv7TpFDMbqeSxLZqJYcxIWpD6eiM0m4wYoaQ89oBqfUx/36xn+3LiUmXW3thrBrF7xRCuU600imvgb33tQz/bhdIZqoCRTOWSoU+/OdUh9/HGJUUTly8EPr6bBdu+CN0+c+l+UkFfHJtuXQtLkDR0524YqBVVg0bXTJzTisFnjPnmc0be6I/bOK4veIZOXUi90iWNIXx9JIUor4ZNtyaNzYjq6z+Qij3MkuNG7M72aNWyHqJMmynbp5YNwTpXLZVh4VdjfScrzBhh3Ln5QiPk2bOy4q/QJJKUKukyTLduoU/5wQom+81MCM0iLJtuUulo+cvPqtRE3VdVHueVvMrLxzLCo0r2oX3s5NdlRDPnWgYgln23L40n97A7WLN6F28SZc3bgpETluvJBk2U6dq6cQfdO87ZC2WPviGrAE6zj+KGcbKtPhOWtbuhWcnzxqsHLoXlIofJZS9/eqUPjMy1/bo22HemGR3M2ajlJu7NJG/JkR5lgs1wBwgT/PcZOWnFVJlu3UKX4g/8Obf/yRjZvgt36xVcbLLw7pZ5n2OYzZhh1O0+FMfQ2mrnmnRx+3HvgIU9e8E0lqXZ2Uet4WL7T+4SOtaUkKi+ROfuWo3adOu7WXZtt7KH0zzdsPpUbxA8mV7VQq/uKIAL9K/+Dq6T3a7HL91w3pF6nAOU2Hs20523oE+45/imxbLpHCJjiTbctpz9xY2ETktKEo6l3qdmkjPj1z3rUKmd+xLnQndT7+ggskd7ILjGAFLGoXb8LIxd19h83brQVv/4nwCr9Y4bQA9N2XnaMzSm2RrFz4XohRN5f2ScdQF70eDemQBhNWLpAgMPK+w4Lyt7MoorY07OoKM7rHZluRhHAxwTthFUEPu2KdkD5Sp/jDUmrN2w7h2iWvh/LefsjU1/i2fgZeKruLhc9JQuEPXfRJyM7XtJM6xR9WDCwD+NP5ZE00/fbm42RUrxIE7ZxN2BhNK6lT/EmIgU06ydimI3ghzhj1JGSLVEXUvh5Sp/hXvf77WK6bui9KSBV2QQVREHW2yEHiioyd1OmzY6fOxHLdNffaFwZPIqW4y7eUiTNMMerQ32W3+7/RpE5hJRSVClw/IaLjRPSuqW0wEW0hon3G/0GmY41EtJ+IOohomqn9RiJqN449YZSpSw1pi4sPuyh9KTBv3jwAuKGcZbsmhrwxgcZSar7ZZKNyA/0pgFuL2hYD+BUz1wH4lfEcRHQdgFkAxhiv+VGhTimAJwHMR75WaZ3FeyaWuiH9YrlukPwpYRalLxXmzp0LAPuKmstKttO2ZnZBnPxacFX8zPxrAMV7qO8A8DPj8c8AZEztLzDzaWb+APkapDcR0TAAA5i5hZkZwLOm1ySeuFIgND84CQP6WmdstGsX1JkyZQoAnCtqLhvZnjxqcGwz2SClTtOWrC2J+HWZDWXmowBg/B9itNcAOGw6r9NoqzEeF7dbQkTziaiViFpPnOhezKAq4tqHcdfi/R+ZsZZlCP/jtPsmNvHz+yI02XaS66gZ2r9PrAn9VmTG4hKfMfnmDZeCP3RrUatfkh3aLWHmp5i5gZkbqqurux37U4Q5xYf27xN7QqimzR2W01uVb0H8/FoJLNtOch21r337kqmRXs+KS/r4TxVWSpvS4sCv4j9mTHFh/D9utHcCMOc4Hg7giNE+3KLdM1EWMUjC4AiyU1n8/L6IRbbtcreHweMJiVAT+YwPv4r/NQAPGI8fAPCqqX0WEfUlopHIL3T91pgynyKiiUbEw/2m13jCaoCEQVIGRxKq9ZQZsch2pr4Gq2aORc3AKhDyM4BVM/XPNqsqe6UuQs0OcWX6RyWccx2AFgCjiaiTiP4GwGoAU4loH4CpxnMw8x4ALwL4PYA3ATzCzAVn9MMAnkZ+UewAgDf8dNhqgHjFrcTd0P59EjM4gkRdSF4TZ2bPng0A1yIhsh0UlfWoVTPHRdCTaJAstP4hTniC64aGBm5tbXU8Z2m2XdnnN3nUYNzdMAIL1u+0PE4APrDI0x8ntYs3+X6tVaEZ4XOIaAczN0R93WK5Lq64BuRdPV8e8WeOhUnMHFw93bJ6VYG6If0SVaQniFwXsKqpIbjLdUlshFuRGascfdP84CRk6mssXTlD+/dJnNIPyrrth91PEmLHruLawQ+7lGS7sOej+cFJlvs/Jo8anCilD0jqhjhJZQUuK1ZkxqLhqsG2ljzQfUNUUkuiWWFXA1iF8wmf0Ql5nCquFUqNfrFxE85Z/Jy9qXvN3LTUXV52+xjH8SqER0lY/AXsLHkgnYXICwSp9VuRnuwBZY1dDQVz+/5V0zG0f59ux4f274P9q9I5S02L4VWKpNLiL665a65cnyZLXpUVmbG+45a/cJlMp9OAauW3JIQY66SCSGalMZA6i9+q5m7jxnYJ7bIhrmymgjc+tolpt2svFWZPuNL9JEE7qVP8dotgpR7aFSRhm5B87PZrlPo+Dok4i4fUKX6nRbBSpvnBSZY5e4TSwG7nbtqyZ/pBjJroSZ3iL1fLCPBfrEPcYMnHbuduqa1XWREk6EJk2x+pW9xdNG205UaXcrCMrhhYhZyPmc2SV9rLQoGknVIMTAibps0d8p35IHUWfzlbRrWX289qnApmf3rGPYWzIKQRP4aQkEKLHyhfy6jlffut+8tnyGYYoTxZmm2XRWKPpM7iL2ecys653QjnrG3R3BtBSAbNkpvfM6L4ywTVRF+CkDTccvrI9i/viOJPEXZlJ6MuRykIurFboxpYVYllt4+JuDelj2iMFLFq5rgeP1gvfJ5jXeKhhbSyfMYYVBZtVKnsRVg+Y4yrG9OtvobQE9+Kn4hGE9FO098nRLSAiJYTUc7UfpvpNY1EtJ+IOohomp6PUD5k6muw5t7x3SKa1tw7/uLAaH5wEgb0tR4EclNQR2Q7ejL1NWi6+4Zust109w1KQRwr75SFXa/4juph5g4A4wGAiCoA5AC8AuA/A/hHZv6B+Xwiug7ALABjAFwB4P8Q0TWmKkaCAm4RTf8jM9YyuufuBv8ZPssNke148ButV44RfkHR5er5OoADzPwHh3PuAPACM59m5g+QL1N3k6brCwaP/XyPp3bBFZHthCMRa97RpfhnAVhnev4tItpNRD8hokFGWw0AczmoTqOtB0Q0n4haiaj1xIkTmrpYHvzxM+tsjnbtgivaZFvkOhwkYs07gRU/EfUBMAPABqPpSQCjkJ8qHwXww8KpFi+3jMRi5qeYuYGZG6qrq4N2URB8oVu2Ra6FpKDD4v8rAL9j5mMAwMzHmPk8M18AsBafT3k7AZiTbw8HcETD9QUTTmFxgmdEthOCRO7oRYfinw3TVJiIhpmO3QngXePxawBmEVFfIhoJoA7AbzVcXzDhFBYneEZkOyE4Re5IxJp3AuXqIaJLAUwF8Lem5n8govHIT3UPFo4x8x4iehHA7wGcA/CIRD3opxDhYFeaUlBDZDtZFOR30YadOHvh8/Y019KOE+KE17tsaGjg1tbWuLshlChEtIOZG6K+rsi1ECZuci07dwVBEMoMUfyCIAhlhih+QRCEMkMUvyAIQpkhil8QBKHMSHxUDxGdAGCXJ+ULAP49wu7YIf3oTpr6cRUzR76NVuTaE9KP7gSW68QrfieIqDWOUDzph/QjTJLSb+lH6fZDXD2CIAhlhih+QRCEMiPtiv+puDtgIP3ojvQjGEnpt/SjOyXTj1T7+AVBEATvpN3iFwRBEDwiil8QBKHMSKXiJ6JbiaiDiPYT0eKQr3UlEb1NRHuJaA8RfdtoH0xEW4hon/F/kOk1jUbfOohomub+VBBRGxH9Iq5+ENFAInqJiN4zvpdJMfXjO8Zv8i4RrSOiS+L6XXQRlWyLXNv2I3bZjkSumTlVfwAqABwAcDWAPgB2AbguxOsNA/Bl43F/AP8K4DoA/wBgsdG+GMD3jcfXGX3qC2Ck0dcKjf1ZCOB5AL8wnkfeDwA/A/BfjMd9AAyMuh/I17T9AECV8fxFAHPj+l3SJtsi18mU7ajkOnZh9/HFTAKw2fS8EUBjhNd/FfkCHR0AhhltwwB0WPUHwGYAkzRdeziAXwG4xTRAIu0HgAGGYFJRe9T9KBQ4H4x8QaFfAPjLOH4XjbIVm2yXu1wb7xW7bEcl12l09RS+mAKdRlvoEFEtgHoA2wEMZeajAGD8HxJB/x4H8HcATDWIIu/H1QBOAPhfxtT8aSLqF3U/mDkH4AcADiFf+PxjZv5l1P3QTCx9FLm+SOyyHZVcp1Hxk0Vb6DGpRHQZgJcBLGDmT5xOtWgL3D8i+gaA48y8Q/UlYfQDeSvkywCeZOZ6AJ8iP/WMtB+Gj/MO5Ke3VwDoR0T3Rd0PzUTeR5HrbsQu21HJdRoVfyeAK03PhwM4EuYFiagS+cHRzMwbjeZjZBTfNv4fD7l/kwHMIKKDAF4AcAsRPRdDPzoBdDLzduP5S8gPlqj78RcAPmDmE8x8FsBGAH8eQz90EmkfRa57kATZjkSu06j4/wVAHRGNJKI+AGYBeC2sixERAXgGwF5mXmM69BqAB4zHDyDvIy20zyKivkQ0EkAdgN8G7QczNzLzcGauRf4zv8XM98XQj38DcJiIRhtNX0e+yHik/UB+KjyRiC41fqOvA9gbQz90Eplsi1xb9iUJsh2NXOtYFIn6D8BtyEchHACwJORr3Yz81Gk3gJ3G320ALkd+QWqf8X+w6TVLjL51APirEPr0VXy+CBZ5PwCMB9BqfCdZAINi6sdjAN4D8C6A/418ZENsv0uaZFvkOrmyHYVcS8oGQRCEMiONrh5BEAQhAKL4BUEQygxR/IIgCGWGKH5BEIQyQxS/IAhCmSGKXxAEocwQxS8IglBm/H8Ds2qfdcd5JgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit,ax = plt.subplots(1,2)\n",
    "ax[0].scatter(glas_gdf_aea_rgi.index, glas_gdf_aea_rgi['h_li'])\n",
    "ax[1].scatter(glas_gdf_aea_rgi.index, glas_gdf_aea_rgi['srtm_height']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29e3Qd1Zng+/vOkWRsI2zht7ElYyAOyBBiGWOHNOEVOuQ6IUAIj3R3shJimEVPT2561g1JOm7G6fTQM9290usOc8EwfTu5g83LJjy60w0mhkCCH5IbYhljMMaWZdmSH7KRsZF0ztn3j6pdqqpTdU6dl05J2r+1wDqvql1Ve3/729/+HqKUwmAwGAxjg0S1G2AwGAyG4cMIfYPBYBhDGKFvMBgMYwgj9A0Gg2EMYYS+wWAwjCFqqt0AN1OnTlXz5s2rdjMMBoNhRNHW1nZEKTUtyndjJfTnzZtHa2trtZthMBgMIwoR2Rf1u8a8YzAYDGMII/QNBoNhDGGEvsFgMIwhjNA3GAyGMYQR+gaDwTCGKIv3jojsBfqANJBSSi0WkbOBJ4B5wF7ga0qp3nKcz2AwGAzFUU5N/2ql1KVKqcX26/uAl5VSFwAv268Lpm1fLw9u3E3bPjNfGAwGQ6lU0k//RuAq+++fA68A3y/kAG37evn6o5sYSGWoq0nw2F1LaWlqKG8rDQaDYQxRLk1fAS+KSJuIrLDfm6GUOghg/zs96IciskJEWkWk9fDhw57PNu05ykAqQ0bBYCrDpj1Hy9Rcg8FgGJuUS9O/QinVJSLTgZdE5J2oP1RKrQZWAyxevNhT0WXp/CnU1SQYTGWorUmwdP6UMjU33rTt62XTnqMsnT/FrGwMBkNZKYvQV0p12f/2iMgzwBKgW0RmKaUOisgsoKfQ47Y0NfDYXUvHlAA0Ji2DwVBJSjbviMhEEanXfwPXA+3Ac8A37K99A3i2mOO3NDVw79XnjxnBZ0xaBoOhkpRD058BPCMi+nhrlFL/KiJbgSdF5NtAB3BrGc416hmrJi2DwTA8lCz0lVJ7gE8FvH8UuLbU4481xqJJy2AwDB+xSq1ssGhpajDC3mAwVIQRn4bBBG8ZDAZDdEa0pm88XQwGg6EwRrSmbzxdDAaDoTBiKfSjmmy0p0tSMJ4uBoPBEIHYmXcKMdkYTxeDwVAIJto9hkI/yGST6+EYTxeDwRAFswdoETvzjjHZGAyGSmD2AC1ip+kbk43BYKgEJtrdQpRS+b81TCxevFi1trZWuxkGg2GUMlpt+iLS5ipglZPYafoGg8FQKcweYAxt+gaDwWCoHEboG0Y0Jg2HIQqmnwxhzDuMXjvfaMe44Bmi0Lavlzse2eRs4K79ztjuJ2Ne6BvBMXIpNKbDMDZZv62TgVQGgIFUhvXbOsd0PxmR5p2oS7Uo3zO+uyMXE9NhiILfPzE+/orVoWyavogkgVbggFJquYicDTwBzAP2Al9TSpVsUIuqmUf9nvHdHbmYmA5DFG5ZNIenW/czmFbUJoVbFs2pdpOqSjnNO/8J2AmcZb++D3hZKfWAiNxnv/5+qSeJuqSP+j0jOEY2xgXPkI+WpgbWrlhmxrhNWYS+iMwB/g/gp8D37LdvBK6y//458AplEPpRNfNCNHgjOAyG0Y0Z40OUJSJXRJ4G/itQD/xn27xzXCk12fWdXqVU1l0XkRXACoDGxsaWffv25T1fVG+bkeiVMxLbbDAYqsuwRuSKyHKgRynVJiJXFfp7pdRqYDVYaRii/CbqrD3SZnfjSWQwGCpNObx3rgC+LCJ7gceBa0TkfwPdIjILwP63pwznGtVBFsaTyGAwVJqShb5S6gdKqTlKqXnA7cCvlVJ/BDwHfMP+2jeAZ0s9l9aE/+7FXXz90U2jTvAbF0SDwVBpKhmc9QDwpIh8G+gAbi31gGGa8GixgRtPIoPBUGnKKvSVUq9geemglDoKXFvO4/s9chom1FXdBl7ujdeRtg9hMBhGFiMqDYNfE652GL7ZeDUY4o/xiPMSW6G/ZnMHv2o/yA0LZ3Hn5Y3O+35NuJrRtOWYdEyHNBgqh1HMsoml0F+zuYMfPrMdgNfeOwLAnZc3ZgnIatvAS03hYDqkwVBZqm0NiCOxFPq/aj+Y9XrBzPpAAZnLBl5pLbqUSadtXy8/2/Cu6ZAGQwUxubWyiaXQv2HhLEfDB2iedVbojB0m2IdLiy5m41W3rX8wgwISxkXTYKgI1bYGxJFYCv07L2+k4+hHrH5tD0rBP72xl5XLm7Nm7FyCPc7LOt02hRUoccX5U/nudZ+ITfsMhtGEWzEze2gxFfoA9eNrASv39WAqQ++pgawZ+8GNu0O1/wPHT1OTTJBOx29Z519yGoFvMFQes4dmEVuhH2SL85tSgr7jfrA1CeH2JY3cvGhOrB6uWXIaDMNPnFf/w0lshb5bMDZMqHOib90PKUh4urX/dEYxe/L4WD5YE4RlMAwvZlPXIrZCH4YEvN70TCaEVTcuzOm3bx6swWAIwqywLWIt9MFakmkvl1RGsfLZdhbMrA99YObBGgyGMMwKewQI/aXzp5BMCKmMlWo/nVF5q9mbB2swGAzBlCOffkXQefMBVt24kGRCAMub56nW/bFMqzyac/0bDIbRQSw1/SDXqtsum8vazR0oLG0/bjvvHq+hZIKvtszhlph5DRkMBkMsNf0g16pbFs1hXG18C4y42zyQyrB2c0dWoRezEjAYomPGS2WInaYfFlgV9w1a7TWkN511UJk7YMwEhlQPfySmicyMN+UYL+YZB1OOwuhnAL8BxtnHe1op9ZcicjbwBDAP2At8TSmVc8o+NZB23DMTAtdeOIO7P3ee88DivEGrJ6X12zp5qnU/6YzyrEhMYEj18AuQlcubWfXCDjMBx5hSx4tRssIph3mnH7hGKfUp4FLgCyKyFLgPeFkpdQHwsv06Jyf7U46mnFbw63eCa6nHddnX0tTAT2+6mLUrlvG96xd4Opqpf1s9/ALkV+0HTQH6mFPqeAkrrWoog6avlFLASftlrf2fAm4ErrLf/zlWGcXv5zrWmeNqGHC5Z2aU5Z7pX5bHfQbXKxI9OY0E89Roxh+wd8PCWWzde8wE8MWYUseLCdIMRyyZXeJBRJJAG3A+8KBS6vsiclwpNdn1nV6lVNaTE5EVwAqAxsbGlv/65GusfLadjFLUJBOgFKmMoiYh3Lp4Lgp4fEsHGQVJge9dv4B7rz6/5GtwUw5b4EiYnMYSxqY/9hhLz1hE2pRSiyN9txxC33XiycAzwH8EXo8i9N0sXrxYtba2Og/rwPHTjoAHEKylHko5NvNyC9NyCesHN+7m717cVdHJaaxTrsl5rAgGw+ilEKFfVu8dpdRxEXkF+ALQLSKzlFIHRWQWEGygD8BtHlm/rdPjEZNOZ7h9SSOzJ4/POVCLHczl2nA1y8vKUi7vDrMaM4w1yuG9Mw0YtAX+eOA64G+A54BvAA/Y/z5b6LG1XW/dtk6ebusklc4gIjTPnuRJuuanlMFcLmEdZpM0mmV5KMfkbDyqDGORcmj6s4Cf23b9BPCkUuoFEXkDeFJEvg10ALcWc3Ct9S+cPYmVz7aTzihWvbAjZ9K1UgZzOTdc/S6mhUxGZnLITTkmZ7MaM4xFyuG983vg0wHvHwWuLeaYQQKv99QAGaWygp6CKHUwVyoeIOpk5C8Ec+viubErBFNtyjE5G48qw1gklhG5QdpwIYI8roM56jV4UjqkFWs2d7BuW6exOftw10N2vy70GOaeGsYSsRP6YdpwoYI8joM56jU0TKgjIYKyVzZRVjdjEbMRazAUTqyEfk9fPw0T6kK14TgK8kLJdw1t+3pZ9cIO0hlFQiCREDK+lA5jkSCTn9mINRgKJ1ZCv/vDj1n1wg5WLm+mvesEMgznrNSGaSHHdX9XCzIdPXHVgul8PJjmhoWzxqxAK4fJL9ex42YGNBgqSayEPlgaW3vXCdZv62QglamoLbtS5oFCvXT8ycC0IEsmhFffPUwqnWHr3mM5PZZGM4Wa/KIKcmMeMoxFYpdPv7YmgcCwJEuKkpSpmORuhSR78n+399QAj921lO9dv4BbF88llTZJo3Il32ppauDeq8/3CPyvP7qJv3txV1Y9Az8mKZdhLBIrTX/GWWfw2F1LAVi3rbPi/tP5zAPFaoKFmB2CvuuOSHbfh4YJdZ4EbmOFQjbxC7HzGz99w1ikrLl3SkXn3oHhs7XmOk8p+XOKten7v6s/a5hQZ3LAR0BP1FqQ57tPYffe2PoNI4mq5d4pJ8PlqZPrPKVogoX4kEe51h1dJ4ynSgTK4dprbP2G0Uxshb4mSOMaLi2slCCvcicEq0lIVgnJsUbU516qwmBcQQ2jmdgJfffABrIEZ9B7lRb8xRy/kJQL7ut1CzXPMdKK6y6azqVzJ49Jk8Nwat/G1j/6MOa6IWIl9HWNXD2wr7xgmpNWWQvOruOns96L20Ns29dL1/HT1CSEdEaRTAhdx0/Ttq83NAFbTUJAhFR6SKgtnT+FmoQwkLYic1999zD3uGoGjyWGU/uOaxoPQ3EYc52XWLlsnuxPDeWcSWV4eWe3E6SECH2nB3mqdb/zXjIZroVVq46u7mBrt3SACNdcOANEWLulI8uF0K/JBwm1WxfPdYLU0umx61Y43DWG/a6ghpGLcc31EitN/8xxNaQSwmBaIeBUzAJIZxSPvv4BaftNAb7aEpx5spiZvVzLP3cHS6czfDyYzvK118fXgmxgMEPGdQz3ZHbzojnD4r4ad4z2bSgWY67zEiuhD4AIoEgkhCSWBqxlf8rORyMK6moT3LJoTuAhCjUFlHP5V0gRbi3IfrbhXX67+wgZlT2ZGWE3RKkbtMauOzYxY8hLrIR+76kBzrDzzmQyituXNHK4r58X3+52vqMUJBPCyuXNZXO1LKe9OKiDLZhZH9rhWpoa+O51n/BMDP7JbDQkmqs2xq47cokSx2KEeXTKUS5xLvALYCaQAVYrpf5BRM4GngDmAXuBrymlchrYj300wCz772Qywc2L5rBpz1E27Ox2TD0KUErRe2og9DiFzuzlXv75hXQ+oW00kcpj3DBHJrkm66gTuZnwvZRD008Bf66U2iYi9UCbiLwEfBN4WSn1gIjcB9wHfD/KAQW46hPTnEhUt91bsDT9oJQJ+vs7uk6ggFsiVpuKg9B1TwxGeyk/xq47Msk1WUedyM2E76Uc5RIPAgftv/tEZCdwDnAjcJX9tZ8DrxBV6Au8squHDTu7ncyT7V0neKp1P6m0AhF2HerL8ufXrpyap1v3s3bFssiCPw5C12gllSEOE7uhcHJN1lEncjPheymrTV9E5mHVy90MzLAnBJRSB0VkeshvVgArAOpmDuW10Ru4OvPkOZPHk85Y76XSGVY+205GKepqLDOQOwe9+xiFzurVFrqjVSuJw+rF7I2MPHJN1lEncjPheymb0BeRM4F1wHeVUh+KRCuBopRaDawGGDfrAkduJxNWuUD3zKxnaxFxJoDBVAaxP/O7PtYms81A+aiW0HWbp0abVlLtidRQXoZ7As81WUedyM2EP0RZhL6I1GIJ/MeUUuvtt7tFZJat5c8CevIfByf4ZuXyZnpPDXg6lp6tdcZJLRhvXjTH2fTVNv2evn6m148r+FqqsRQMKqTiv/aRzGhdvYxFzAQ+8imH944A/wvYqZT6e9dHzwHfAB6w/30237HmTz2Te69f4Ag7rVHA0EytO1iQG6T+d83mDp54tp10RhVceasaS8GgQipRUziPBIxNdfRgJvCRTzk0/SuAPwa2i8ib9ns/xBL2T4rIt4EO4NZ8B5pQl3SEXT5XrVx+uyufbSdl+3gODBbeMYtdCha77B3tQnG4JtI47BuMdkZ7Xx0LlMN753UIrWF+bbHHDdMo8i0vN+05SsZVGCYR4N5ZCUpZ9o6FjaZK21TXbO7wbO4bs0NlqERfrWb69EKJa7sKIVYRuW7CNAr3ZDAwmOFnG97lu9d9IjufTSpDQoRVNy6siolmuFYXYxV/SmrP6s6YHSpKOftqkLIEw5s+PSqjZT8jtkJfaxTrt3V6XDH9Scp+u/sIW/cecx5AtbRm/yQ1VuvZDgf+wXfzojlOIj6AhAzP6s5QOmEZMOO4bzBa9jNiK/Q167Z1MpDKsN61IetPUuY3/xQq8P2/advX60w2UaJ69e+1142pZ1tZ/INPgHG1liKQSAzf6s5QOmEr+jjuG4yW/YzYCX23AA6bWYOSlC2dP6XolMp+d8n7n7cENuSP6g06Z7EaQdu+XtZt60SwUiobwRWMf/C5XXbDJvvRYIsdjYStzOO4xzVa9t5iJfT9lbNWLm8OnVmDHsCDG3cXLGz9AvpX7QcZTA2FeOWL6l23rTOrklcxGkHbvl7uWP0GA2nLTPFUWydrv5M9aRnhFT74CpmYx+q9iyNBewRx3eOKa7sKIVZC/2R/ihqfv/pjdy11tF8//gdQjLANyn+/+YNjjqafK6q3bV8vT7cN7TkkE0LDhDqPqSeqcN605yiD6SG7dNCkZYTXEO6kW+7XQYwWW6zBUA5iJfTdlbPcmTTX23b9fIFWxSy/gn6zYGZ9JJv+pj1HSaWtyUGAqxZMz2vLD9PUl86fQm1SHE0/aNJyC6/+Qet+jFXhVcgEOFpssQZDOYiV0P94MI0oK2c+du6eQrW0QpZfbgHsjoDVx/BHBPvxC5Np9eNytjWXoGppamDtimU5bfr+QulPt3VGTh892iikX4wWW6zBUA5iJfS7jp9mpu16l0oXbx+PQr58N1E0Sb8wAXLWs80nqPJNWC1NVqH0NZs7UAwVSh+LG5eF9ovRYIutFnHuS3FuW1yJldB3++NrX+tStLRcHcIT5JXKZEVzRtUk/cLE31Z3G8oxgUUplD4WbP/l0N6NwMhPnPtSnNsWZ2Il9EUgAUhCuPqTQ+n3o2pp/ijNXB3CLYD9qZpLWWG4Nxh3HerLsvGXKqiiCLuxsnFZivZuBEY04tyX4ty2OBMroT9/6pncdHkjT7Xu5+Wd3bz23uHIgzEoSjOfKSUsVXOxHjgwlAMmnbE2ozNKedpw79Xnl9wx8wm70bZx6a41UK6U00ZgRGPp/CnUJK2+lEzGqy+Ntn4+XMRK6E+oSzJ78nhSGRUp0jZXIJcurBKlQyyYWZ81ARSqAa7Z3METWzvYfuCEU8RdC35BDWunHE0bl3oy17EQCaEsmrkRGLlxT7TYCQwzGSsyHnK7yJZyvkL662jq58NJrIQ+BA/GKEmZ/IFc+aI03cesSSb4asscbrG/X6gGuGZzBz98ZnvW+0k7JUA1CqIUa/qIm51bPw+931MuzdwIjHDcYyMh1mpVAakMPLa5I3KNiqh9qdQMtebZFUbshD5YvvFuH/mgSFvwJmXSgVxRozT9G7lrN3ewfltnzijgMH7VfjDrvRpb4N95eWPxN6IAyiGs42jn9ifYSwgkkwkOHD9N277ekgV/ta8vjrjHBkplfZ4ri6l7hRB1xWxMbcNLrIS+Pw3DLYvmANGTMhUyiPUxtdnAXYS9UA3whoWzeO29I87r6y+awd2fOw8gb6bNcgjrcuWSL2bwlWtlEHYc/97LK7t6eHlnN49v6fAk4SvX+Qze8ZZMiJWy2iX7w7KYBq0QovSlfKY286zKS7lq5P4jsBzoUUottN87G3gCmAfsBb6mlOrNdRx3GgZ/rvxyJ2XSx1y/rZOnWveTzqisyWPN5g5+tuFdblg4K6fGrj/7VftB57tRtOZyaNZZlcJK0JQKtXO77e3JElY2+e6DO1hu5bPt6GwVxV5rvvONdSHjHm9dx0+zZnOH81lCCM1i6l8hJCLuZ+UytcVx9TnSKZem/0/A/wB+4XrvPuBlpdQDInKf/fr7uQ5y5rgaMnly5bspZnnuH9AtTQ0e2z9Y2nnf6UEe+s0eAEeLzyf43Z/7teb12zqzOnU5lrWb9hwtWy75Qu3cm/YcdVZKqYxi5bPtLJhZX9Q1RLkP5brWXOczQsbCPdGu29bpSVvtHwduk45baSjE+y1sLBvTT/kpi9BXSv1GROb53r4RuMr+++fAK+QR+hPqkjxs58p//b0jjr09lwdPIeQb0G6/ej+/aj9YkBbrWSInEzzVup9URnmif/2DpBgBtnT+lLLmki/UROYs/4GMyp2RNNdxotyHYq/V329ync89kRVTX3k04a8TERRwCNkOFeV0XDBeVuWnkjb9GUqpgwBKqYMiMj3fDzRn1CY93hp9pwfzpk3IR9u+Xn624d3ADeEgO6Q/q+cNC2dFbT7g1ZoPHD/N41s6HLOV2/7uvg7IvweQ6zzlGGiFTKwtTQ2sunGh53qKGZRRr6GYaw2b6MOO0zChbqjv2a/HIkGmOy3w/ftufocKdx6rUon6zMe6Sa4Qqr6RKyIrgBUAs+bOczqam0df/4C+/pSnNm6UjUu/J4Hb11trDbnskN9cNo8dBz/Ma9MPw71EXm+nThDfBpceJFHNCn4tKyhhXLEUY9q48/JGFsysL3nARV1h+L+Xb7CHmQfCztd7aoCEWMpGQmBH14kxWfbSb7r78S+3O8/ZfT8VXoeKSpQJzdc3jEmuMCop9LtFZJat5c8CeoK+pJRaDawGaFxwsXL7ZGt0igR32oR8ngF+TwJ9DAGm14/j4jmTgezlY7mXp5A7+tctvPPZLj2xBQkBEVLp8nX0Yu2n1XJ9jDLY/c83n1DKZZYbS8Jk6fwpJARn0zytYNXzO7jtskbP/bxl0VB8S7XKhObrt3FeBVSjbZUU+s8B3wAesP99Nt8PMrZg9qOAhbMnBXYuRHh1Vw8Hjp/2pBn2aPBYkbFa8B/6sJ9Db3fz6q4e1q5YltNzoJwJvfQxgjTjKMLJ07nTlqOpO19QqZ0mrvbTsOcQZZIKmnCjZk51m+Xc5sC4CpBy0tLUwLUXzuDFt7ud997qPMGOrnbu+uy51I+vddKVaBPYxl09WVXkhuMe5eq3cV4FhAWdVrp/lctlcy3Wpu1UEekE/hJL2D8pIt8GOoBb8x3n8Ml+ZgVI/QTWstvNlRdMY8Pb3aSVYsveXrbs7fXUs/V3hG8um8fq1/Z4Yk10KcSgfDi6fOFgWlGblJx1coPIlzu/GOHk959GhHS6fAK6lP2BNZs7PC6r5SLsPrbt6+XA8dPUJBN574G+34WW01w4e5In70zDhLrYCpBy07av1zGFuhymSGUUj77+AatuXOgxmfpxF0Eq5tyF9EF/v4WhvbE4e//427ZuW6dTMKqS/atc3jt3hHx0bTmOX1MzNODCOhl469n6O8KmPUezggsV4Rt167Z1OlWsBtKq4CpVxXQ23e4fPbM9UGMK6tzlWom4j1XI/sAD/7KTxzbvo68/DURzby2EoPsIeMxcty9pjFRIPspKxp+eI5Ox95eUor3rREHPtJpmhVLO3bavl9se/h3aiU2wMuC6c0qt/s37OcfiVQumF90n73hkk/OMgupEB+HeP4taZ7va+PujwLBMUFXfyM2HAF9tmUPvqQGC7P1u/PVs/Rr1uNoEH7s2iROSvYJwnzfX63wUaypp29fLU637h+ru+jIb+q+p2E7hHlzJpJAoYn/ggX/Z6cQyuCnUvTUXQffRPRGkM4rZk8cHrtSCUnIUkpba7bqbzqiCkvhV06xQ6rkfevV9XJeOAi5ramBbx3HHRLr36Kmcx5hWP66otmtNF6z7v74IZcvtctvedYKbF83xVKOLi40/SInLVyujHMRK6E87cxwJsdJ9KCyzTl3tUDoGdw4WPwJcOndy6LH1DfZH4IZt7N28aA5PtVkPIJEQmmdPKvh6/DmEorBpz1HH711PeO6OGjW9cL6O7R5cqbRORGENNHckdK7j/uuOQ4HnLtS9Nd81BG2u1yQTlnDGcul1fz/MPBZlsC+dP1SSUpOASEn83FTTrFDKudv29bJhZ3fW+5Mm1PHE3ctY9fwO3uo8kfMYtUnhZnvMFopfqVMUtmrxu9zqsa7Tra/Z3MGP7dTndQWsJCqFX4kbjiSAsRL6iYSlT2tb4hXnT/UIn8fuWsq6bZ08sbWDtE/yK2DLXmtZettlwct9fYNvjuBt0NLUwP1fanZy4696YUfkaNMgX+Z839cPWmu2OgDprHE1/PCZ7Tzd1kkqnXGW2HU59hmiaHr+wZVMCMr2iHrtvSO88f4Rnrj7M57frdncwY9/uZ2Msgb2Fy+elaXx3XPl/LJo+fmuIW13AKVwVhv/9Mbe0LwvEF5Uxy9U3CUpg/phlD7gfo4iMqz+/qVsyK/b1hmUY40TpwZoaWqg+ZxJOYX+p+ZMYuWXmovWqG9ZNIcnt3aQykBNwtpXKWTV4na5FSyFRptJ12/r5PEtHZ40Hg+/+j6r/2RxUW0t54ohrF53JUhU9OgFcua4GupqEiTtnOl+bbOlqYFzJo/3dErx2V1SGUs4ff3RTbTt66VtXy8PbtxN275ez3Huvfp8x2TkFw6a3lMDTlrZoM/DCLNDB6GF29+9uIuvP7oJgJXLm0nYka4P/WYPazZ3uDyRLAbSilXP7/BcVyHnXzh7EvYcS01S+MmNC5k96Qzn81TGWua72/njX24nba/CBtKKruOnuefK+Zw9oZY5k8/gr2+6mPu+eGGke5SPXNewac9R0j7B9Ms3D3hMPqBt0eJ4mQQdz3//tbCvzdEPo9DS1MDK5c2IWO25/7n2wGdVCfSq9nvXLyjYtBNmxtRj6ZZFc7KERsL+ryYh3HZZo8e2ru9rIdeeSFj27UQiwQ7XPkr/4FA+/zAaJtSRELFicZLiPMfamgQ9ff1Z/ebld3oceVFIW0u5vkoeKwqxEvoT6pJ5O6vWYpICZ9QmuPsP5lueLC60kF63rZM7HtnE3/7bLu54JPtmejpIgEaU7/MwdBstB5vcWl6QMNrRdcIx8eTirc4TgdflvkfJhNBlpyHWuDV2NweOf+x5/dvdR/jRM9sdLcQ/YLbs7eVfdxzi5ECarhMfs+qF4EmoGNzX4L/3S+dPIemTTj19/dQkrXuusDQ9HVi06oUdTsoL//Hc979/MMNf/HI7j2/pAKW4fUljTqEZpFC4ae864Zkk123rzPubcqEVm0InqzCzjALHTPRXN11MUqwJoi4prLhyPgnbJXrls9ylhVIAACAASURBVO2s2dyRc9L23wP3ax3EqMdwT1+/M74Vlrkm7N617etl1Qs7nLHzrSvOZe13LHmycnkzr757OOs3mYziV+0HIytpmiiKVdRnXYiSGHYOqRs/MepvYmXegWB3Rr+XSZCt3L+hmEgIWz44FroppDuIrm61crl3SapNP/rzby6b5zyMKF44K5dHMw0F+ef/w4Z3Q4/td6EL2uzSmt66bZ083dbJ2i1W4YuVy5vZ0XWCtZs7PPsiqbTV8f3TzKmBNI9t7uCptk6+9Zl5ge1xm3f0ROte8hayBPZ/N8y+2dLUYAlkV/ZHpaz9j/3HTvH6e0c815IrZbbbDJNhKH38YDp4g9jd1nxmh93dfVmvh3NzN9+994+rdds6OdLXTwKy9s1qXC6Y7qyyAjxuB7DBUOK9VTcuDDQx+dOAr1zezP3PW8GKNUnLxKifnQJe3tntWdkPpBUPvfo+l86dHBi3oaP5M8qK5P9880zuvfp8Hty4m5TfJoylFN2wcBZb9x7L69W1blunsyG8dH7uMpKFbKaX4vShz1HbMPsTkX5EDIW+G32jn2rdTyqtSCQArI6hbeVt+3p55PUPsn6bySh295z0vOcWBO5d/nTGcsdz++ZrjVHZx3rk9Q8sDw6BGz81mwtm1OcUZH7TkF8YaoJcS3Np+QKcP20iuw9/FHhd7uNu2nPU2Qf4eDDDj365PdBeC7B5z1Ek5FiDqQw7Dn4Y2ibdLhGcYKaapHDNguls3NVDKq2oSQq3Lp4buqkdNkjC7u/Ni+bwVOt+Z8M1mRQEaxN5856jDKaHIrC1z3jQ8fT9/9mGdz01EfzuvP6N9K7jpz1ePv7Nb20OcdOfygzb5m6U9NFut9cMekM/mAzw8KvvM61+HM2zJ7HqhR0eTzjPd5XyTLLavLbrUJ8nDfjHgxn+8bcfOIrZYMD5g4bCS2938/LO7qzrWjo/PAGgFqxuN1N3oaNcqUS0XNB97am2Tu7/UvOQhhAwqNZv64wUqJbLWSEf3gBUiWy1ia3Q153S3bGsiXrIy2TTHivfdzqgZ/j7T0IsO7bmve4+j0bxVOt+Dvf1Ow82rSzblzYj6HMoBb98swuwzEthD8sfSKUnLu17DF4fe/dv3b+7asF0Xt7Z7VyPApbMn0LHsVNO4FjYRrHfEyVM4AMebxU/iQR05HHRu2xeA637ep1BmkorTzTnYFqxZnN44ZNcHidh7pdrVyxztNNX3j3M2i0dzkrI/WwRYdehvtBB3dLUkFUIB4bceR/4l52sfm2Pc22CXcHLPllGwevvedOAb9pz1COwkra9e1d3dgqOSuB3PfULHX90dz5jovt56k3/IHQNY3e07srn2i2lTbKF+Ps+xSwqQf2kpSk8AWBLUwPfXDaPh22LQFLgrs+eS++pAdZs7gjNIqqfpXtCGkxl+FX7QVK2+2o6o7L6ay63a02prrUNE+qG7mcBPuWxFfq6U4Yhdi71dXk2djQZBfc/b5lZXtpxyBHcmsG0Ynvncd9J4LYljRzp6/cIME2/nfhNm4B0eLrfPPHW/uPO7wdSGb6/7vd0HDsV6BcfZNbQbmaZjKIm6c134o9A9K8ibl0812MGKYZ0BvYdyy30O46dCtTK/ISlK9Zul8KQzV2v9LTnUpD7pb4PG3Z2k1HBmuFAylrlaD/7oMHlj9dIJoQ39x/naw/9ji17vRq7wlYqMoqLz5nE7ztPoLD6gw7iWzo/Ow10uZLT5ULfl77Tg869yKjsIES9X4WyVmEpe/KKQjqjrNWBT/B/5VJrBdwwoc4x2cDQBBx0/IinDCRImC6YWc/XLpvr8csH67488toe53xpBQ//Zui1YMXxrFzenOXRt3T+FGqTQ8pTIiE0zzoryySk733X8dOBbtd+1kVcDYTRe2pgaHVewI2MrdD321r9XPPJoYi/J1xuWG7G28FY+iPttrV2S7AQPPRhf9Z7Qo4ALhlKBqc9bbSm4zZP+Iumu81OQQVW9H96k6bv9KBj60ylM7y045AzuYA7NTQ0z57EsvlTnM+bbS+dqAO6WILuXRBB6YrXbO7w7Ml8c9k8gKwI7DD3y5XLm7N86/3o2I+BgPtNQJsyGcVLARO9m7SChedMYufBDxmwteWnWvc7Jqywms3DYdJJiDgCwZ/GRO9nZeyssl+8eFaWEuTH3YcSAosarZgYPSEmBf542TxamhpY8YvWnApb2PGt8RT8GQw9Q8gWpm4FQa+Sm2dPcj4PckRQvr8HbA1e97l+u3rfDQtn8dXFc9nd3Udbx3EydioKt5IHOKbhZAJPehD3StxtJny6rXNoNVBE2gqtWFiTq4p8w2Mr9N2Dpu/0IG/sOUr7gROOj/g9dg3aXYf6mDnpjCzPE4CJZ9QwkBpwHraIZbvOJwC1fbommfD49fq55JxJtHd96LHB+5edbft6OdwXLhCTyYQTd+DO8aMHsD/UPaOsTWutmbjzmWeU5dHzVucJS2NOiuPJEid2dHn9vP2F5Xcc/JD68bVZEdjJZIK+04Pc+1ibY/YbSFlRl5bvbv4LzSgCC9q8td+7yot6y3r6+rl07mRH+A2mlbOxXkkBH4TbZJNR1uoTlV2ucP22Tuf+KaX4TYBXi59Zk87gotmTeNleUflXP2kF3/lFK+dNm1iUZ5JSeFZN2lrh7r8JhkxLyYQ1qelz+U3B/kpu2uMrh16AAEf6+j2mwdffO8Jr7x1x9oYyLiXvkdc/4JpPTqfr+GmPaTiVgc9/chrT68fR09fP+m2d7DrUR3vXCWfVqmNJ9HlvXTw3p9NAvv3A//j3XeHeHz5iK/RhyEtGpwuoSQqL5k5mIJVh16E+dh3qy9Ki3Rzps7Qbsf+XUXg2QMP45Mx6ln9qNgd89UH97D58kk83TmabPfvrzqrdNN2CO4izJ9QyflwNB3pPA94cP3oAh/VRrZkorJWFf1NNEc1WWw3Wbu6gefYkp3rYx4Npz+c6otc/WaVSmSwvrYyyBuqgT7OsS4Zr/vp97aJZyqS44e1uRxPVuA9XiaCfsM8bJtT5TLvKk5dIa8NPbHV5PQG9pwbztunA8Y9pnj0p577QsY8GOPZR8Ko4H7U1Cc+eh4hkOTTUJIX7v7zQEZ7aK+2WRXMCx1g6ozwb7D/5ysVOrImfhAACOw95Pa7cE0DGduTQ9yDtWg3644UEPI4G+j3nla92R5CrbC7zpkYrFn86cDq/YLOJldDv6et3Zm69BHpia4dnd3+rrWG81bmd86efGem4hdq8dh3q469umsKuQ32B7muak/1ptu7tJSlwyZxJnDt1Ii/8/qDjpnnlBdNyJqU6dmoQfAPuiL0qyB7A2WSUtTl9qj8VuDyPo8AH636GTdZfuXQ2C2bW82drtwX+Logg81suU49GkXtzOwrK/i9pKxXujXW/uSVK4fhCvG78exyrXtjhEWjpDLQfOOGkH9D7T/5LjnoLdnSd8Ai9cvKtz8zz7Hk0TKjzePoANE6ZyIKZ9fSeGnC80rRPv9trR6Owoss3f3CMtd9Z6hz/oVffZ8Pb3R5TkV5lhCFYKWGWzDub3/g2/MG6J1qoJxPC1PpxWX3Q/aomKXzrinOdIk1+pwV/4ScIrxXw0KvvUztl7idz3V83sRL63R9+zB2r3wAR52Hm4vhp72AXLNe9XK5nfkTgkzPqPTN8BviLZ7bzwdGPfNNzMGnbrNLe9aHj5dM/mOHld3oKFrwvvd3Nil+0snFXT86lqOaZf++kde/wRHoOB8++2cVzb3aFCvggWvf2DvsEZy33dYWt7E38tn1WeU49aDNK8cNntvPKrh7u/tx5oVp/vrw5YZ+7TTZu3uo84diaS71HXcc/do6REPjyp2bz3FtdZTEfbtjZzX1fvNBjEus4+pFnZbe75yR3PLKJ+7/kzZx5y6I5LJw9yfHaqUkIjVMmOntn/liWjQHjctn8KVYwXUjHSySES86ZxOu7swW+xtkkzih+a5uEwm5NWin+8Xd7SaUzbLYDMptnTxraAA+YXLXdX08M73X3Ocqe1I4bucFZVueM1ouO9nmFvqW5FdYDlYL3AtzG/Mu8KLhdR5XvdeT2QKCnUBhbKyDwh2PjNwytORf6m0oxoTbBqQBhqoB5U8/kg8MnSWcU//TGXo/W7bcxa158u9tyLw1J9JUvUCfo87Z9vaHOCRBt1RMFz8anggnjavJuoEdl9+GPWPGLVs+EWD++Nut7A6kMO7pOBG6S33bZXCdo8+FX3/c4TOgWWq602Rq4v9aGn3RGZe1j5CKft1s6A2k7bfeA7c7srlQW1KkVVmqUV989XPBGuZvYCf1Cuk/Qd8Nm6lxESXkwljC3Y4izJ9ZxKsBJALxeWP2DGf78yTdZceV59J4aCN3HgdzuebkikcM+X/GL1oo9szAFQGGZIt2uiaU24cW3u9mws5v5087kW1eca01wAXszT7Xu5+ZFc5zEZH6T18LZk3hl11B11pqksHD2JB7cuNtJxxHkIFFNFLk3mcFSiPN5lEUhVkJ/fG2y2k0wGBwEOPRhsMD3CzmFlZLih89s554r5+c+ruBZpgd5ZfjttjoITTPVzle/ZnNHQSvDQgkThgm7Df760n/74jsc+yj/xnCu8+3uOencx68unsuRvn72HD7pOGEM2skGb7usMSs62h04BdZzumbBdE+qh6sWTEeAX7/TTQkK84il4kJfRL4A/AOQBB5VSj0Q9l2/F4fBUE3OmzaRPUeCnSJybWhueKeHpikTQguNLJhRDwzFISQT4Zu8OjlekBYYFp9STnRtaTc6FkVHuOtAKIATp4sX+H5Wv7bHaUPT2RM80daWa/J2EmKlU3D7xU+ZWOe0a1ytlZ3A7Qyy4e1uxtUmuOuz83mydb/lUDGGqKjQF5Ek8CDweaAT2Coizyml3g76vrEqVJdxSaG/0lJkhCDAdRfOsPI6BUj3XOaAwyc/5sSpVOjn7xzq46FXh8oNap9yICsdwMpn20MFe6Uf1fnTz2TJuWdnuS1fcf5Ummed5Ul30Dx7Er9qP0imjJqzE1GcVqGu1hll2dtvW2KtCHYc/NDjyfbNZfOyJm6FFRn+qJ1Pq1R0Pq5n3+waETKs0pr+EmC3UmoPgIg8DtwIBAr9AtJHGCqAEfhDKHCS7BVKLoGvj/3yzm6PgNCCXwtRbbcvVijpsZRIWHNW0GHOnlAbquXWJoW/ueUSdvkcGnRWSn/itB89s72sAi9XnIUfsYsNBcXUvLHnqBW850N7VCksU9VZE2o5HnIv8gV1XdbUwM9u/zSnBtIVNbWFIUDmdF/+CDubSufTPwfY73rdab/nICIrRKRVRFrPqkln5cY3GEolIVZFp1p/Ev48lEMLDCPo0Dr/jc7ho3O+FMo9V87n6f/wGf7zHy7gO5+dH6hMLZnXwKnBdKii9e0rzqWlqSE7BkIpdnSdyLo35b5TuQS+v82ZjOKNkBz0A6lMqHNHTUJIirXR+4mAmB/B+szfFF2nQvPm/uO07evl7s+dR9InUS+cWc/1F83gzssb+eubLubrlzfmVG515G8hT/2nN11M6kR35ARbldb0g9ruuYVKqdXAaoDFixer6y+by1q7epHBUA7q7GjPlc+1V7spkVBYKaoP2ykBCvWM0ZkkP988k1sf+l3gBKPTPIcdV6fS1platVafVvBud19gMJSbcnjzhB130vgajp8eWk1lFIyrydZfm86ewK7uYNdrhZX6QAFPt3WydW+v02adwuTWxXPZ/MExj5eWAD/5ysW023UptGu23pC//bJG3u3uYyCV4bbLGrP2aR7cuDtvgJtOaLf8klmR4iDCcoOFUWlNvxOY63o9BwjN7HRqwNrIra1JhM50l81rYMm83OHs91w5n3X/4TOc0zC+sNYaRiVfaJ5J76kB0lUyX/lD9KOQUVagXjFBVQorP9OfP/lmoMBIiBWMlEuY6FQYLU0NLL/EW+h+695eq7ZEyG8TPk24nCjwCHx3m/zsy5H5NSHW5vM5k8c7gaDKfv+OyxtZu2IZNy+aw96j3v2Au+0a0LcsmsO42qFKbA0T6vj6o5tYu6WD7QdOsPJLzYEb8zrOwn93pp5Z51wfWIL/ghn1PHXPZzh/WnjcVcL2BCuESgv9rcAFInKuiNQBtwPPhX15z5GTTqm6S+ZMyvq8JiHc9Ok5PHnPZ7jnyvmcdUbwQuX9Ix/R0tTA5z4xLeuzqF2xiFW1IaY8+1YXb+4/XpTwLQfu2sPDSWdveIBQX3+24Dx/2kQ+NWcSf33TxY7AatvXy/NvBaf4CJKnAtyxpJFVNy7kjNpw5a3cFDoxrviD+U4iNne/yCjrGpyaCK5Z4/qLZjg1oHW8hC7tmq/etkb/7o7LG6lxCZkjJ33ZBWxh3tLUwOUhQj0h8FdfubjgpH4VNe8opVIi8qfAv2G5bP6jUmpH+PeHduObz5nEru4+O/kSKKysdLr04H1fvJD7vnihlXPk+R285cqb8et3epwizk+37neKjdz/5YX0nhqg7/QgD+eIwDt/+pnOJpZ7g0qH3Vci98hIJpmASWeEbwpWG6WsxGhRmHZmHb2nBwtK5ZGP8XXVCYc5K+SZZBQc7uu3Upfbrow1Cfibr34qS4A89Or7BXsJzZ483slz87MN7/Lb3UeqHvykmXnWOP7s2k/kzIGkm+qPfr7bzuyrCSt+lK9Ajv6dQGitCz0pgV0lrq2TATtl9Hd8KT8KpeK9USn1L8C/RPmuCM5yyV0opOv4adbaZfiCquWs/FIzX3v4jaHqVraN7d6rz2ftimVZATAPbtztsTmePbGO46cGnIRZf3PLJc6D2dF1gjV6j0FZbnwvv91dUG6Y0cz1F83gqgXT+Ytfhmc7jQNRZc6RkwMeDaxUEgJ7DhdXHapUcmXPVMD9X2rmx3YStkQie9Hftq+XX7/Tk/W+AA0hnj9ugdfS1MB3r/sEW/cec+piaB//LzTPDHVxDNoP+Mqls3l5Zzd9/cXH8tQkhAe/3pKVy8itxCUFJ2FevuhoN4V8V3Pzojk8sXV/1t7IVy6d7awo9LHXfqewY+ciVhG586eeyVeWDO1uu4uJrNvWGTqLtjQ18BNdJi2jqKv1djz/TdJFjbWWc+wja6DfsXiup9oOWA/Gfe57Pnce0+vHeWZo/8aMncY8sENrIdnedYItvk2ifAiWC16hqSbsrLFZE1VCLNODwtIK3znUl1M4Xjiznnd7TjqVwn5iBxQ9uHF37FY/OiHaNZ+cwSu7eiLbxhXBtVqLQZu1K6HlCvD5i2Y4FcOCyHXajbt6mFY/zslVNZDKsOr5Haz8UrOn8Egm4OAKKwGZn+svmpGVTM4tDHV9YS24lpw7xcn86W+ve1P1zsut9ND/vN1bd6FQ7vrsuYGywF/hLFd0dC4KrZ+gyzu603sLcIEdvFfKsXMRK6EPVoGHAbuQuLv6lL/IMuC5CYWUomtpauCrLXM8fr2ptGL25PGRcp0Anolg5fJmJ8e3jgo8d8rErKRt91w53zODt+3r5bbVb0Q2JYjAbZdZy9Ld3X0cOPExvR8NOBvggb/Bikr85rJ5PPLaHmep7hba7vas29YZOhl9uqmBv7rp4qx7oQfOx4MZBMs9UoGTC2XZ/Ck8+vqeYQ15V1gT76fmTubuz53Hum2dPGnXKc5FubxOBGuJ/k9v7M2ZXtvPhLoEpwczeSfR86ZN5O7Pncf+Y6dyJgdMYAtPgen145wKZ6m0Ynd3nydh2ludJ7jjkU1OMjh/oXE3fht0Qqx7HZZPKOh995jtOz3Io69/4GTJRMQZSzfbK/5SXWiDErgVo6GXEz3+ghTWShEroX+yP0VNQNpYf47psFzjhcyGtyyaw5NbOzyCyF8yL+y4YR3FXbd23bZOz2C8bF6DR+Dr4zyxYllWfu8wdH1cfU/uWP1GqD+zXm0kE+JUh9LfFKyMhH67pntlpRNY6Y0tnSc+6B63NDWwcnmzE1z0TnefU2u0fzBDe9eHtDQ2MMmuETCtfpxTRMWpimanpU4IzJ92JnsOn/RoP+JL/DVvygQ+c/7UQPfepGvl9db+4zRMqONwXz9nnVHj5IXRx9QrsoRY2uvUiXWRSz+CZQv3T2Z6g+3Oyxv5fPNMHn71/chBO589f1qk7+4+/BFfe+h3ee3tSizBn1HZJS237u3l8xfN8Bawd427lqYG7vrsuVmFa9xofb8mpPh3PpxCSavfIJVRJAXu//JCFsysZ/22To99PWpGzwTw5Uu9EbJ1OezsQbmOhnMSGI7ayW5iJfTPHFdDxt4MSSaELrtylRb0MDRIiykk7KalqYHbLmt0zDQJKczfNUz4ud9zbyLfd8OF/kM4v3nkTxZz60O/87idnTP5DMbXJnn/8EfOMtddE3TTnqOhA0AXrba0XeUsqf05yHNdW9iSPIzeUwNklHIm7Ce2djiphXVa2jpXOUg3D27czfYD1ka8UnDu1Il0HDs19MyxwtyftwvUAHQdP81Z42qc3DCOf3VNgi8utOzFGWVlbgwTogtm1LOoqYHm2ZPY0XWCp1r305OjtKWfZEK4YPqZnsn9U3MmeUwkLU0NfGruZF6KMKlrt71kRBOe//HrokKelMIKwtaBCjh+asAT/eo3nwZpx/5jiD5Rkazb1umcP62gvesEC2bWW+/bufAfu2spVy2YnndCdE+4f7xsHuu2dWYVSc9FviI2laKc5pt8xEroT6hL8vBdS50SYWu3dJCQoUGtSUh25ywGv73ef7xSZvyWpobATeQwPjGj3iP0r1ownZsXzeHrj24KFNRBNT8TWD7GzbMnseqFHZ7JEyhoGVtoJ3RPKslkIqsOLli2cl2+DnDa0jChztHiFZb3VUvjZE/+8ud/f5BrPjndWRGl0srJnaJzn1wwo56+04NWbvQ87VVYNRPe7TnJbZdZ7+mI2KikMyrLtLLMdW/dK9Ra21MmIbC4qYGtAYVfjpwc4MW3u4t2c2yYUMu/dxRWX6E/lWHtimWhwtHa/8pdmMgdoFSM4PJfr5BdLGbdtk5PuuSEWE4V4K0/cfuSoYCoYgTppj1HHXPcwGBpimVciZXQhyH/WF0ODeWtSylYCZ903ctSzxUmCMsx4xfS6bRblhbwevCFta+lyVvzMylWpKDu8FpTctcSfeyupU4O8nLjbuuB46eteAsfCqvQ9OY9R0HEqfvprw+azijauz70vJfJKKbXj2NcbcKpoaqVAaXghd8fZPklFJz0Kp1RrN3cQW1Nwoo8TSsyeAVRIcfTkazu8oTWpqd1lIzKX/jGfz73PoMA4wMKuwiW4pCrKPmFM+sZTGc8yctuu6wxZz9taWrga4vn5qwVrVdYxSphzbMnZb1eMLPeszIVvHUvMgp+895hVi5v5jfvHY60go1Cw4Q6515nCDf5jmRiJ/Qh2z/2m8vmWZs89kZHOQS+JqzD5ytbV27C3LJyDcggW6B7dXLO5PGeWqLDcQ0tTQ2s2dzh0ZjnTD6DM1ymKss7RjlmuiMBJhX/5rSeCPWmnqPR2+dJZ1RgneAoKCCdznD7kkY6jp1y/MqTArfZ3mT+Itdh3LBwlpMd00ldUGKeGsVQDV4FgZW8Fs9rcFauYRvH7/acRClFTVJonnVWYJqAIKz9r2zXQrAE/mcvKE0J6z014ExsCfu1X+EBsq5tMJWh99RAlmJUygq999SAk8K5UJPvSCGWQj9Iw/1888xh3VzJV7auEhSzHHX/xr86Wbm8edivAbyDGODA8Y+prUlQmxTHHKMQlFLU1iSYVj8up9eMYPmUu+/N1x/d5JlYwn7rL2zvf61d8mtdK46te495NMcW2+7vL9St23betIlMHFfjCNEHN+6O5GnivuZ8JSpVjmvUx9Ib6g9ufI8DAdW+nDalFWeNr2XBzGzXwCC0a+GPXe6V+pzjyqCEae8vfz91OxZs2nM02wPN3jzONQaCVui5JoVqjPvhJpZCH4I9ZobTtlZtV65i8K9OgrSgQilGa9KDWGtlWqP//EUzmFo/jqdtM5b2LNKmqDANVXwal75ON0FCM+EzC543bSIdvacd2/qKP5gfqEw8Zu8ruU08vacGsgS5Nmv4I1ndLqxu5kw+g06XML5sXgP9qQwzzjqD+VMn5vSSyTeHHPtowIpOf2FHYG1e9+ZwBvjt7iNs3XssstnyzssbPUnGEsAVJWr4mqhm1oSIx8zldmzQ5Fuh55sURuK4L5TYCv0oVNq1argnmlIJ0lLc2tKDG3cXdK+K3dfQA2fdtk6esj2YFPDKrh5uXTyXVDrj8Szyewtt3NXjcWH1uwPq6/RPElrwZ7ljAnW1CS6fP4X3tW1aWZ4pYc/YHy8S5LOuINBrRWvc/qCbKxdM52k7nD4h8O/7j5POKHZ19zHNLn9YLOdOOzNwMgTL3dYfcFaMye+WRXNY73J8GG4zK8raHxFUqP0+n6YexWw70sZ9oYxYoR8kkIBRPUPnI0xLKVZ4l7Kv4R447hS0iuAcJe7v954asIqMqGCNzj2puAPi9N5PKqMcWZyQoY3/Xa6I41ybdEHXvXT+FMd7yC1Ww7xWek8NeMP7E1acw8LZk5zN94wtiAdTmYJcRf0kE3CPnRemribhpDzQBEUYF7P5Wg0t2C/EdcxJWJBmvjaOBfNNPkas0PcPzPXbOh2/3uH0r40bQVpKscK7lAGiV2ELZ0/y2GvdOZXCBEeUmAJ9ne5jbdpzlIxP865JDmmkm/YcjbRJ5z+/Tps7kMpkpbIQsWIG2vb1Bpp4/OH967Z1ZvnX19YkmJ5nX0NTf0YNfR97M2Re+8kZHtPUpj1HeWv/8VCf9qRYro1RfdfdxMHMGsVEk8sbabSbb/IxYoW+f2AqcASbrjwU5YEOd/RdNdDRjINpK2dOVOFd7AAJ2lD2B3jlOpZbkxdg16G+0Db4B7hb0/UHDTVMqCMhlmjNF6Hpvm6PicGHQjwuse7rC7p3WT7pgmdfYyCVCc3bBPBRfyprcpjqMg25G8ymJQAADqpJREFUPaiChH7C59o7EvA/41I960a7+SYfI1boB7l0PW271Smsaji35NFkqhV9VxVsYVdoUvliBkjQhvK9V59f8L6CtqtrO/242tzPSPcJdzpfbX4BWPXCDidZ3MrlzXknHv9kErTRrDd3g4SP36tk056j1I/zDjllR6DeeXmjpz/vOtTHE1s7+P2BEx4zUUZZpQ7b9vU6WWGDVkJhq5gVfzB/RAn8IIyJpjRGrNCH7IF5qx1Eov2u82kAw+2LXwzlWInoYLeo96VUggZloROsX7vWHkD52t7SNJTO131+fTz3BnJU9MZskMtm0nYRyiV8/B4ofp5q3c9COxeRe/P9zssbWbO5gx/9crtH8J8/o57v33BhJBOZe2NXyJ9WYSRgTDSlMaKFvp98aRX8xF1jKNdKZLivM2hQPrhxd0ETrFtoZZTlgRO17WFCQZt+RKTgSEudW8jPNZ+czqVzJ4duLEKAB4rPvTSVVk70btK2/2ttXP+rN3/rciS+89+Dtd9ZysOvvs/L7/SgVLBJa6SaN8e6iaYURJWQKElEbgXuBy4EliilWl2f/QD4Nla+pz9TSv1bvuMtXrxYtba25vtaTgrtxHHu9A9u3M3fvbjLiQ793vULik6jUMx1lvPe6AlMTzxRJjB37pqwpG+FtFGnRsjYArCQSTQs8+jaFcsAck7O/msPSnPtjgGoSQhP3L0s6xjFPouw344p8+YoR0TalFKLo3y3VE2/HbgZeNjXgIuw6uE2A7OBDSLyCaVU8WVvIlKoBhBnjaGcGnqh11lugVDMkjxfmwttoz8TaL7Vhl9YhmUezbeKCbr2zzfP5KFX36fnw4/tegMfOKajjMp2Ay2ln4b9diSYNw3lpyShr5TaCSDZdsobgceVUv3AByKyG1gCvFHK+cYa1bRdVkIglDrB+oVwUBt120v10Q6bUKK4mAYdN+i3r75rJQrbeaiP5ZfM4rm3ulAqd+73chJ386ahMlTKpn8OsMn1utN+LwsRWQGsAGhsHNleBZWgWiuRuAmEICGcy5++1BD7Qia9YiZn7ZkElqvx878/6Cl6E2ba0W0r1dRVbLsNI5+8Ql9ENgAzAz76kVLq2bCfBbwXuHmglFoNrAbLpp+vPYbhIW4CIUgI33v1+aH+9DoHeyFZS91UetILc/30exa59yF0GUGdkto9qZWSMqPaz9YwvOQV+kqp64o4bicw1/V6DlBc3ltD1YiTQAgTwkH+9DqZ29NtnYECMgqFTHrFCNxbFs1xKqslE5BIJJx0Evra/Cma/Smp3asPY583RKVS5p3ngDUi8vdYG7kXAFsqdC7DGCCKEHZ/p+v4adZu6ShJCEad9IoRuC1N3spq+jjua9u0x1sMPCFWOuFUOtvtNG7mOEN8KUnoi8hNwP8NTAP+WUTeVEr9oVJqh4g8CbwNpIB7h8NzB+LtgmmITtBzLGTlUT+uxkm5UA4hWIkc7P7rCTquLrOYTAg/uXEhgOPTv+qFHSyYWe8cJ07mOEN8KdV75xngmZDPfgr8tJTjF4rxOx4dRMndk+t3Ol2CELwxWmp7hjUHu1LWdQhOlTRd9N6/qoiTOc4QX0ZVRK6xa44O3M9xYDDDj59tJ5OxNPa13wmfyN2pFoCiUi7ka89w5mDftOeolSaaoRxCxoxjKJVRJfTNgBgduJ8jDHm2DNgptMOEq86iqWxNOCGlFewOak+l+5XbjBRWFMeYcQylUFIahnJTjTQMhniin6M/L/ydlzfy1zddHPh9bdpJJoS7Pnsu9eNry9YPhqNfmcJAhmIZzjQMscPYNUceuTZt2/b18squHgbTKjSNMHhNO0op6sfX5sxTVEwgU6X7VVgsgunPhnIy6oS+YWQRZZPU7doYJgDLkWIh6HvVLA1ozJOGSmCEvqGqFLJJmqsIS7lTLFTDE8zY6w3DgRH6hqoSVbuNIoTLmWIh6sRQbgFtzJOGSmOEvqGqRNVuy+mOG+WcS+dPoSZpp3RIBhcfueORoRz5uVxJDYY4YYS+oepE0W7Lbe+OpFFrz7YADzd/lsxcrqQGQ5wwQj9GGHfTcIbb3h0UGOU+p38aiI/js8GQGyP0Y4JJIZGfcti7o06s+VYW7iyZuVxJDYa4YYR+TDApJCpPIRNrvpVFVFdSgyFuGKEfE4yPduUpdGLNt7IwnjaGkYgR+jHB+GhXHjOxGgyjMPeOwZCL4d4sN5vzhuFgTOfeMRhyMZwmGbM5b4gjiVJ+LCL/XUTeEZHfi8gzIjLZ9dkPRGS3iOwSkT8svakGw8giaA/BYKg2JQl94CVgoVLqEuBd4AcAInIRcDvQDHwB+J8ikizxXAbDiELvISTLlNffYCgHpZZLfNH1chPwVfvvG4HHlVL9wAcishtYArxRyvkMhpGE2Zw3xJFy2vS/BTxh/30O1iSg6bTfy0JEVgArABobG8vYHIOh+hi3TkPcyCv0RWQDMDPgox8ppZ61v/MjIAU8pn8W8P1ANyGl1GpgNVjeOxHabDAYDIYiySv0lVLX5fpcRL4BLAeuVUP+n53AXNfX5gBdxTbSYDAYDOWhVO+dLwDfB76slDrl+ug54HYRGSci5wIXAFtKOZfBYDAYSqdUm/7/AMYBL4kIwCal1D1KqR0i8iTwNpbZ516lVLrEcxkMWZjgJ4OhMEr13gmtPK2U+inw01KObzDkwgQ/GQyFU6qfvsFQNUzwk8FQOEboG0YsJvjJYCgck3vHMGIxwU8GQ+EYoW8Y0ZjgJ4OhMIx5x2AwGMYQRugbDAbDGMIIfYPBYBhDGKFvMBgMYwgj9A1jgrZ9vTy4cTdt+3qr3RSDoaoY7x3DqMdE7hoMQxhN3zCqCNLoTeSuwTCE0fQNo4YwjV5H7g6mMiZy1zDmMULfMGoI0uh18JaJ3DUYLIzQN4wacmn0JnLXYLAwQt8wajAaffUx9Q3ijxH6hlGF0eirh/GSGhmUWi7xJyLyexF5U0ReFJHZrs9+ICK7RWSXiPxh6U01GAxxxnhJjQxKddn870qpS5RSlwIvACsBROQi4HagGfgC8D9FJFniuQwGQ4wx9Q1GBqWWS/zQ9XIioOy/bwQeV0r1Ax+IyG5gCfBGKeczGAzxxeypjAxKtumLyE+BPwFOAFfbb58DbHJ9rdN+L+j3K4AVAI2NjaU2x2AwVBGzpxJ/8pp3RGSDiLQH/HcjgFLqR0qpucBjwJ/qnwUcSgW8h1JqtVJqsVJq8bRp04q9DoPBYDBEIK+mr5S6LuKx1gD/DPwllmY/1/XZHKCr4NYZDAaDoayU6r1zgevll4F37L+fA24XkXEici5wAbCllHMZDAaDoXRKtek/ICILgAywD7gHQCm1Q0SeBN4GUsC9Sql0iecyGAwGQ4mU6r1zS47Pfgr8tJTjGwwGg6G8mNTKBoPBMIYQpQKdaqqCiPQBu6rdjghMBY5UuxF5MG0sD6aN5WEktBFGRjuD2tiklIrk/hi33Du7lFKLq92IfIhIa9zbadpYHkwby8NIaCOMjHaW2kZj3jEYDIYxhBH6BoPBMIaIm9BfXe0GRGQktNO0sTyYNpaHkdBGGBntLKmNsdrINRgMBkNliZumbzAYDIYKYoS+wWAwjCFiI/RF5At2la3dInJfFdvxjyLSIyLtrvfOFpGXROQ9+98G12fDXiFMROaKyEYR2SkiO0TkP8WtnSJyhohsEZG37Db+l7i10XXepIj8u4i8EMc2isheEdluV6hrjWMb7fNOFpGnReQdu28ui1M7RWSBfQ/1fx+KyHfj1Eb7nP+nPWbaRWStPZbK10alVNX/A5LA+8B8oA54C7ioSm25ElgEtLve+2/Affbf9wF/Y/99kd3WccC59jUkh6GNs4BF9t/1wLt2W2LTTqz02mfaf9cCm4GlcWqjq63fw8oS+0JMn/deYKrvvVi10T73z4G77L/rgMlxbKd9/iRwCGiKUxux6o58AIy3Xz8JfLOcbRyWGxzhQpcB/+Z6/QPgB1Vszzy8Qn8XMMv+exZWEFlWO4F/A5ZVob3PAp+PazuBCcA24PK4tREr7ffLwDUMCf24tXEv2UI/bm08yxZWEud2us53PfDbuLURS+jvB87GCp59wW5r2doYF/OOvlBNaKWtKjFDKXUQwP53uv1+1dstIvOAT2Np0rFqp202eRPoAV5SSsWujcDPgP8LK1OsJm5tVMCLItImVqW5OLZxPnAY+H9tU9mjIjIxhu3U3A6stf+OTRuVUgeAvwU6gIPACaXUi+VsY1yEfuRKWzGjqu0WkTOBdcB3lbdecdZXA96reDuVUmml1KVY2vQSEVmY4+vD3kYRWQ70KKXaov4k4L3heN5XKKUWATcA94rIlTm+W6021mCZRf8fpdSngY+wzBBhVG3siEgdVv2Pp/J9NeC9SvfJBqwa4+cCs4GJIvJHuX4S8F7ONsZF6Me90la3iMwCsP/tsd+vWrtFpBZL4D+mlFof13YCKKWOA68AX4hZG68Aviwie4HHgWtE5H/HrI0opbrsf3uAZ4AlcWujfd5OezUH8DTWJBC3doI1eW5TSnXbr+PUxuuAD5RSh5VSg8B64DPlbGNchP5W4AIROdeehW/Hqr4VF54DvmH//Q0sG7p+f9grhImIAP8L2KmU+vs4tlNEponIZPvv8Vid+Z04tVEp9QOl1Byl1DysPvdrpdQfxamNIjJRROr131j23fY4tRFAKXUI2C9WUSWAa7GKKMWqnTZ3MGTa0W2JSxs7gKUiMsEe59cCO8vaxuHaOImwgfFFLC+U94EfVbEda7FsaYNYs+i3gSlYm33v2f+e7fr+j+w27wJuGKY2fhZrCfd74E37vy/GqZ3AJcC/221sB1ba78emjb72XsXQRm5s2ohlK3/L/m+HHhtxaqPrvJcCrfYz/yXQELd2YjkVHAUmud6LWxv/C5aC1A78f1ieOWVro0nDYDAYDGOIuJh3DAaDwTAMGKFvMBgMYwgj9A0Gg2EMYYS+wWAwjCGM0DcYDIYxhBH6BoPBMIYwQt9gMBjGEP8/BJA9e5ybEE0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "glas_gdf_aea_rgi['diff'].plot(marker='.', linestyle='none');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.559596869346697"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glas_gdf_aea_rgi['diff'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### geoid offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 581206.91326115, 5394760.93117372,  589576.38352026,\n",
       "       5411418.92006629])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glas_gdf_aea_rgi.total_bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fd86cc822b0>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALIAAAD8CAYAAADT2P50AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAXu0lEQVR4nO2dfbhlVV3HP18GEOUl0BmVAEUSLS0DvU4Uj6Xgy1iGZmloFonFk1mPZD01xOMLpk9ovmSSL1NOoCGIoL2KCCgYycvMKOiMgDMC2QAJo6hYvM093/5Y68zsuXPvPfvss+8565zz+/CsZ/bZe6919r1879q/tdZv/X6yTRCMO3uM+gGCoA1CyMFEEEIOJoIQcjARhJCDiSCEHEwEYy1kSask3Sxpi6TVo36eYHRoXOeRJS0DvgE8D9gKrANeYfvrI32wYCSMc4+8Ethi+xbbDwLnAy8e8TMFI2LPUT/AABwC/Hfl81bgZ+beJOkU4BSAfffRM378qU9fsgfasGHDNtsrmtZ/waqn+Tvbfljzu269xPaqpt81aYyzkDXPud3sJNtrgDUAz3jyPl6/fv3SPZD0X4PU/862H3Lt+r+ode+eetXyQb5r0hhnIW8FDqt8PhS4Y7EK2v8nl/SBBsWYTmd21I8xloyzkNcBR0p6AnA7cCLwytE+0oDYdDoPjPopxpKxFbLt7ZL+ALgEWAastb1pxI81EMZ0vH3UjzGWjK2QAWx/BvjMqJ+jPYxDyI0YayFPHiHkpoSQS8LGnRByE0LIpRE9ciNCyEXRwbP3jfohxpIQckHYYSM3JYRcFIawkRsRQi4Jh5CbEkIujTAtGhFCLgi5g7bfP+rHGEtCyEURpkVTQshFYRSmRSPGeYfI5GGgM1uv9EDSYZK+IOlGSZskvT6ff6SkSyVtzv8eVKlzWt7/eLOkF1TOP0PS1/K1v5E0ny/4SAkhF4VRZ3utUoPtwB/b/gngGOB1kp4CrAYut30kcHn+TL52IvBUYBXwgbwvEuCDpF02R+ZS3M6UEHJRuLUe2fadtr+cj+8FbiRtD3sxcE6+7RzgJfn4xcD5th+wfSuwBVgp6WDgANtXO+1U/milTjGEjVwSNtpe27F+uaTqvq01eVvXbkg6HDgauBZ4jO0709f5TkmPzrcdAlxTqbY1n3soH889XxQh5JKwa/W2mW22Z3rdJGk/4CLgVNs/WMS8XWgPZK29kaMmhFwYanHPnqS9SCI+1/an8ulvSzo498YHA3fl8wvtgdyaj+eeL4qwkYuiPRs5zyx8BLjR9nsql/4FOCkfnwT8c+X8iZIelvdBHglcl82QeyUdk9v8rUqdYogeuSBkt9kjHwv8JvA1Sdfnc38OnAlcIOk1wLeAlwHY3iTpAuDrpBmP19nuPsxrgbOBhwMX51IUIeSSsNH2B1tqylcxv30LcPwCdd4OvH2e8+uBomMphJBLI+JaNCKEXBRGnc6oH2IsCSGXRHeJOuibEHJR9DWPHFQIIReGHKZFE0LIJWHD9odG/RRjSQi5JGyIwV4jQsiF0eYS9TQRQi6K6JGb0tPXQtJaSXdJ2lg519oug7y2/4l8/trsctitc1L+js2STqqcf0K+d3Ouu/fgv4oCMEnIdUqwC3Wchs5m9x0Bbe4yeA1wj+0nAu8F3pHbeiTwZlJekJXAmyt/MO8A3pu//57cxgTgEHJDegrZ9heB78453eYug2pbFwLH5976BcCltr9r+x7gUmBVvnZcvnfu9481stH2h2qVYFea2sht7jLYkZ0pR6H/PvAo5s/adEi+9j3vDJJW5I6FxkRv24i2B3tNdhn0W6evHQvV9GSPe9zjFrqtDLo2ctA3TR3rv53NBVrYZbCjjqQ9gR8hmTILtbUNODDfO7et3bC9xvaM7ZkVKxqnwBsSYSM3pamQ29xlUG3r14DPZzv6EuD5kg7Kg7znA5fka1/I9879/vHGQMf1SrALPU0LSecBzybt2t1Kmkloc5fBR4CPSdpC6olPzG19V9JfkNKQAbzVdnfQ+WfA+ZLeBnwltzEBGLZHpKEm9BSy7VcscKmVXQa27yf/IcxzbS2wdp7zt5Cm5CaLbo8c9E2s7JVGeL81IoRcFGH/NiWEXBJhWjQmhFwaIeRGRICWgrDB212r1EHSquy8tUXS6iV+/JESQi4JA52apQfZWetvgRcCTwFekZ26JpIQcmm0JGTS9OQW27fYfhA4n+SgNZGEkEvDNUsOK1spp8xpaSGnq4kkBnslYXCndlaDXmFlxyIcbFuEkEujvfWQhZyuJpIQckkYvL01a28dcGR23rqd5MPyyrYaL40QclEI6psWi5I3KfwByYtwGbDW9qZWGi+QEHJpuL3MX7Y/A3ymtQYLJoRcEv0N9oIKIeTS6MSMaBNCyCVh4dkQchNCyKURPXIjQsgF4bCRGxNCLor2pt+mjRByYbjF6bdpIoRcEiZs5IaEkItCdGLWohEh5JKIHrkxIeTCiFmLZoSQC8LEYK8pIeSSsMK0aEgIuTDCtABJ+wL3V+IG9iSEXBIWnl3W+74JQ9IeJMf/3wCeCTwAPEzS3SQ31DW2Ny/WRrzHCsMd1SoTxheAHwNOAx5r+zDbjwaeRcqAcKakVy3WQPTIBTHFg73n2t4tMUoOI3wRcJGkvRZroE56ssMkfUHSjZI2SXp9Ph8pytrG09kjzydiSSf0uqdKHdNiO/DHtn8COAZ4XY5YEynKWkfYe9Qqk4Skl84pvwqs6X6u00ad9GR32v5yPr4XuJEU6CNSlC0Bnt2jVpkwLgBOBl4E/HL+d9/KcU/6spHzK/9o4FrGJEXZuGV1mjSzoSY/S0rnsQ74kG1LerbtV9dtoPaftqT9SIb3qbZ/sNit85wbWYqyccrq5Ck1LWyvA54H7A18XtJK+oyKVOs3kkeMFwHn2v5UPj02KcrGiWEM9iT9laSbJH1V0qclHVi51tpAva+f2+7Yfh9pLvlP+q1fZ9ZCpKxJN9p+T+VSpChrG6fptzplQC4FftL204BvkOZvWx2oN8X2HbZfbvuIfurV6ZGPBX4TOE7S9bn8IsmmeZ6kzaTXwpn5QTaRjPevA59l9xRlf08aAH6TXVOUPSqnKHsDeQYkzyN2U5StY/cUZW/IdR7FhKQoG4aQbX+uMr64hp1vyjYH6rWRdISktZLeJmk/SX8naaOkT9bt4eukJ7uK+W1SiBRlrWJEp/4S9XJJ6yuf19he0+BrTwY+kY/bHKhv6+MZzgbOI5mU1wD/ALyV9AZeS5qhWpRY2SuJFsPKSroMeOw8l063/c/5ntNJ6wTndqvN/1SNBur9sL/tD+Zn+n3b787nP5Lj1/UkhFwYbS1R237uYtfzKumLgOOzuQCDDdS3zhmo90NH0pNy3UdImrG9XtITSQEYezJZ8zgTwDBsZEmrSGOME2z/X+VSmwP1fvhT4F/ZaXuflsc+XwLeWKeB6JFLwkPzozgLeBhwaR6XXWP799rMJd4Pti8Hnlw5dZWk5aTZkFo+ySHkgjDQ6Sy9P3KeKlvoWmsD9bpIeibw37b/J3/+LeBXgf+S9JbKTNWChGlRGB2rVpkwPgw8CCDp50lTuR8Fvg/UmomJHrkkhmdalMaySq/766SpxK4f8vV1GogeuSC6jvVDWNkrjWUVd4Pjgc9XrtXqbKNHLowJFGkdzgOulLQNuA/4D4A8/fb9Og2EkAtjGoVs++2SLgcOBj5Xmb7bA/jDOm2EkAvCFrNTuIsawPaOpXFJh5AWQu4HbqlTP4RcGNPYI0s6DdjL9lvzqauB75H8k88mO6QtRgi5MKZRyKR56GdVPn/H9tHZhfRKQshjhpnEOeJa2P7fysf35XOzkh5ep34IuSDSVqepFPJ+kvbqbvm3fTak3SfAAXUaiHnkwpjSeeQLgQ9LekT3hFL8tw+xc6f8ooSQC2O2s0etMmG8kbTn81uSNkjaANwGfJvwfhs/7Okc7GUPt9WSzgC6Dk1bbN9Xt40QclFMpENQTyTdAFxF8j/+T9u39dtGCLkwprFHJoUA+DnSJuY3Z/v4S91i+9peDYSQC2MahWx7I7CR7LKZnepPBE4F3kWN7U4h5IKwmcSBXE/ywsfRpF75WFKs5NtJoSOurtNGCLkoptNGBn5ACo75t8DqHFOjL0LIBWHSosgU8jukQIa/A7xa0jpST3y17dvrNBBCLowptZHPI/kkkxdFVpJMjL+UtLftx/dqI4RcGFNqWnRX8n6GnXbyM0kRjP6zTv3pG1kUTb3l6bZ6bUl/Isl5lqB7bujROCV9BfgWKb7FMuDdwOG2j7YdkYbGjWHOWkg6jDRv+63KuWo0zh8FLpP0pLzy1o3GeQ0pZdgqUmyLHdE4JZ1Iisb5630+zknA1xoEdtlB9MiF0UG1Sgu8l9QDVsUzkmictr86iIihXnzkfSRdJ+kGpaxOZ+TzkdWpZYa1i1opY9Lttm+Yc2mhVBeHUDMaJ2mz6KMGesAG1OmRHwCOs/3TwFGkZDTHEFmdloB6wVnygHC5pPWVcsouLUmX5RjDc8uLgdOBN837ALszjGicA1MnPrKBH+aPe+Vi0ivl2fn8OcAVpMB4O15PwK05JthKSbeRX08Akrqvp4tznbfkti4Ezsq99Y6sTrlON6vT+aSYua+sfP9bSH8oY00fve2iYWUXisYp6aeAJwA35BfiocCXlfJ2jCoaZ/fZlgG/BBxORZtzMiXMS90cIstyxJe7SMLaLasTUM3q1NbrqZWsTt1e6+67767z444MA52apfF32F+z/Wjbh9s+nPS7e3qOuzaqaJxd/hX4bdL/3/0rpSe1Zi3yqPUopaQpn5a0WzC7CsVldSI7o8zMzAz9ldcXI/a1GFU0zgqH5rwmfdPX9Jvt70m6gmTbflvSwU459trK6jT39bSVneZLt84VVLI65V55crI6DXmJOvfK1c9Dj8ZZ4WJJz7f9uX4r1pm1WJF7YpR2tD4XuInI6tQ67m+wN4lcQ3rj3yfpB5LulbRYTscd1OmRDwbOyYb4HsAFtv9N0tXABZJeQ5pUfxm0+3qy/V1J3axOsHtWp/MlvQ34ChOS1alTtvGz1Lyb5DzU9+JInVmLr5J8Reee/w6R1al1ptT7rctmYGOTwWIsURdEWqKeaiHfCVwh6WLS+gVQb/othFwYLS0/jyu35rJ3LrUJIRdEd4l6WrF9RtO64TRUFNM5ayFpTV5xnO/avpJOlvQbi7URPXJhTOmkxQeAN2YxbwTuBvYhTd0eQBrsn7tw9RByUZjp3CFi+3rg5ZL2A2ZIU773ATfavrlOGyHkwpidQiF3sf1D0spt34SQC8JTHB95UELIhTGlNvLAhJALI3rkXak4hi1KTL8VxDD8kUtE0lWV44/NuXxdnTaiRy6KiYxGX4d9K8dPnXOt1i8khFwQZmpnLRYbGtQaNoSQC2NK3TgPlPQrJFP3QEkvzedF2mTRkxByYUynjrkSOKFy/MuVa1+s00AIuSCmdR7Z9qsXuibpMXXaiFmLwpjGWYu5SPqR7Ch0GfDlOnWiRy6IKR7sdfeDnkCKVfJ0UhiAl1DTtIgeuTBSirLeZZKQdC7wDdLm4rNIAVrusX2F7VovoBByYQwriKGkP8yx+TZJemfl/NDDypL2cd5DSr9wU96s3Nefawi5IJIbZ70yCJKeQwpT9jTbTyVlTmo1bl8/5LiCLyf5Hl8m6T+A/SU9tm4bIeTCGJJp8VrgzByfD9vd4DojCSubn+Em22+y/WTgj4CPAddJ+lKd+iHkoqhnVmTTYtFonD14EvCsbApcKemZ+XwRYWVtr7f9BuDxwGl16sSsRUHYMFu/t100Gmeeuprv1Xw66f/7QcAxpFwdF0g6ghGFlZX0/h51ruzVRgi5MNpaol4orCyApNcCn8pmwnWSOsByRhdWdn3l+AxSTOy+CCEXxpBm1v6JFF/6CklPIsWQ2EaKwfdxSe8h5RDpxu2bzXHYjgGuJcXte39uqxu372oahpW13bWxkXRq9XNdQsgFMcTNp2uBtZI2Ag8CJ2XxjTqsLDT8Ww4hF8YwFjtsPwi8aoFrowwr25gQckGkJepRP8XwkXQvO3viR1RCyYqU/eOAXm3Unn7L6Re+Iunf8ufI6rQETKPTkO39bR+Qy56V4/3riBj6m0d+PWkJsUtkdWqbmqt6U+p8vyh1k+EcSsq28/eV09UVnXPYdaWnrdWhHVmdbN8DdLM6iTTqvnCe7x9b3EcJdqVuj/zXpCyZ1bdaZHVaAqJHbkadHCIvAu6yvaFmm8VldbI9Y3tmxYoV891SFNPoxtkGdXrkY4ETlBI+ng8cJ+kfyVmdANReVie0e1an+drakdVpnrbGFgPbXa8Eu9JTyLZPs31oTmN1Imnl5lVEVqclIWzkZgwyj3wmkdWpVbr+yEH/9Jsw8gpy2M/I6rQEhP3bmFjZK4xJW+wYFiHkggjTojkh5MKYRl+LNgghF0bYyM0IIReECRu5KSHkwuhEl9yIEHJhhIybEUIuCNvMRo/ciBByQcT0W3NCyIURg71mhJALo8+d9EEmhFwQMf3WnIj9Vhi2a5VBkHSUpGskXZ93z6ysXBtFWNmBCSEXRHKsd60yIO8EzrB9FPCm/HlkYWXbIIRcGK7538Bfk2IRQ9qN091dM7KwsoMSNnJh9GEjL5dUDf63xvaamnVPBS6R9C5SZ/Zz+fwhwDWV+7qbeh+i5sZhSd2Nw9vq/yiDE0IuCGM69XvbQcLKHg/8ke2LJL2ctLvmuYworGwbhJBLwu35WvQIK/tRUsAdgE+yM17JqMLKDkzYyIUxJBv5DuAX8vFxwOZ83ObG4aESPXJBGNg+nJnk3wXel3vQ+0mzEa1uHB42IeSiaKW37f0t9lXAMxa4FmFlg8FIK3uxRN2EEHJJCDqKReomhJALI3rkZoSQC8KYWWZ73xjsRgi5MMK0aEYIuSDSyl4IuQkh5MIIITcjhFwUqU8O+ieEXBAmbOSm1E2Gc1veHXB913Uw0pMtBWaWh2qVYFf6cRp6ju2jKq6DkZ6sZbqDvTr/BbsyiPdbpCdbAkLIzagrZAOfk7RB0in5XKQna520JFKnBLtSd7B3rO07JD0auFTSTYvcW1x6MmANwMzMTNHrv8lpKHrbJtTqkW3fkf+9C/g0yV6N9GRLQNrs1LsEu1InYeS+kvbvHpNShG0k0pO1jmPWojF1TIvHAJ/OM2V7Ah+3/VlJ64j0ZC1jOg77twmaplhjMzMzXr9+fe8bGyJpw2I7m3ux17L9fNDDd9uEMS93/++1A33XpBEre0XhmJFoSAi5IFJ85BjINSGEXBI2HcdArgkh5IIIf+TmRICWwrA7tcogSHqZpE2SOpJm5lxbcoevpSCEXBRDW6LeCLwU+GL15BAdvlonTIsBOHz1v+/yee/HPnHeoCf9MGhvW+87fCPAPNFfdzh8Abfmef2Vkm4jO3zlel2Hr4tznbfk+hcCZ811+Mp1LiWJ/7yl+JlCyA2ZK+J26GuHyCBhZReizbCyCzl8LQkh5IIwptOpPWvROKys7YWW84fh8LUkhJALoy2HoMXCyi5Cm2FltwLPnlPnigbPVIsY7JWEhzNrsQhL7vC1VA8ePXJRDGcXtaRfAd4PrAD+XdL1tl8wRIev9n+mcBpqztwB353nnMoDd25unAhG2st7Ljuw1r3bZ7eF01CF6JEH4LYzf2mXz3rHizYM1mKkjGxKCLkoTGfHNsSgH0LIxRE9chNCyKURbpyNCCEXRcR+a0oIuThCyE0IIReFh+I0NIlM1TyypHuBm1tqbjm751t+vO0VTRuU9Nncbh222V7V+7bpYNqEvL6tRYQ22woGJ3wtgokghBxMBNMm5EEdz5eqrWBApspGDiaXaeuRgwklhBxMBFMhZEmrcpyGLZJWq8DkPsGA2J7oAiwDvgkcAewN3ADcDiyfc987gdX5eDXwjnz8lFznYcATclvL8rXrgJ8lbbS8GHhhPv/7wIfy8YnAJ/LxI4Fb8r8H5eODRv07moQyDT3ySmCL7VtsPwicDzxinvtGltynzR92WpkGIc8XX2EZZSX3CQZkGpyG5ttDd6HtkwtK7hMMyDT0yPPFatgMRSX3CQZkGoS8DjhSKeXv3sArgcugnOQ+S/WDTxWjHm0OowC/CHyDNOPwLtIsxA3AJlIIKUg27OWk3vpy4JGV+qfnujeTZyby+RnSH8E3gbPYuVK6D/BJ0sDwOuCISp2T8/ktwKtH/buZlBJL1MFEMA2mRTAFhJCDiSCEHEwEIeRgIgghBxNBCDmYCELIwUTw/0Ro5yv76sPZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots()\n",
    "im = ax.imshow(dem.read(1),cmap='inferno')\n",
    "glas_gdf_aea_rgi.plot(ax=ax)\n",
    "plt.colorbar(im,label='HAE (m WGS84)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "module 'pygeotools.lib.geolib' has no attribute 'get'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-43-78209e892e66>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mgeolib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m: module 'pygeotools.lib.geolib' has no attribute 'get'"
     ]
    }
   ],
   "source": [
    "glas_df['geometry'] = glas_df['geometry'].apply(geolib.get_NED)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "https://github.com/dshean/pygeotools.git"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = gdal.Open('/home/jovyan/data/baker_1_m/baker_2015_dsm_7_wgs84-adj.tif')\n",
    "ds = gdal.Translate('/home/jovyan/data/baker_1_m/baker_2015_dsm_7_wgs84-adj-sub.tif', ds, projWin = [-121.855968, 48.774327, -121.796142, 48.712009])\n",
    "ds = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gdalwarp -r cubic -t_srs EPSG:4326 /home/jovyan/data/baker_1_m/baker_2015_dsm_7_m_utm-adj.tif /home/jovyan/data/baker_1_m/baker_2015_dsm_7_wgs84-adj.tif"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
